[37] | 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 | // Start measuring the time |
---|
| 139 | osmapping::os_GetTime(&tim1); |
---|
| 140 | |
---|
| 141 | // Create the integral image |
---|
| 142 | Image iimage(im, doubleImageSize); |
---|
| 143 | |
---|
| 144 | // Start finding the SURF points |
---|
| 145 | if( bVerbose ) |
---|
| 146 | cout << "Finding SURFs...\n"; |
---|
| 147 | |
---|
| 148 | // These are the interest points |
---|
| 149 | vector< Ipoint > ipts; |
---|
| 150 | ipts.reserve(1000); |
---|
| 151 | |
---|
| 152 | // Extract interest points with Fast-Hessian |
---|
| 153 | FastHessian fh(&iimage, /* pointer to integral image */ |
---|
| 154 | ipts, |
---|
| 155 | thres, /* blob response threshold */ |
---|
| 156 | doubleImageSize, /* double image size flag */ |
---|
| 157 | initLobe * 3 /* 3 times lobe size equals the mask size */, |
---|
| 158 | samplingStep, /* subsample the blob response map */ |
---|
| 159 | octaves /* number of octaves to be analysed */); |
---|
| 160 | |
---|
| 161 | if( bLoadRegions ) { |
---|
| 162 | // Load the interest points from disk |
---|
| 163 | loadIpoints( sRegionFile, ipts, bVerbose ); |
---|
| 164 | } else { |
---|
| 165 | // Extract them and get their pointer |
---|
| 166 | fh.getInterestPoints(); |
---|
| 167 | } |
---|
| 168 | |
---|
| 169 | // Initialise the SURF descriptor |
---|
| 170 | Surf des(&iimage, /* pointer to integral image */ |
---|
| 171 | doubleImageSize, /* double image size flag */ |
---|
| 172 | upright, /* rotation invariance or upright */ |
---|
| 173 | extended, /* use the extended descriptor */ |
---|
| 174 | indexSize /* square size of the descriptor window (default 4x4)*/); |
---|
| 175 | |
---|
| 176 | // Get the length of the descriptor vector resulting from the parameters |
---|
| 177 | VLength = des.getVectLength(); |
---|
| 178 | |
---|
| 179 | // Compute the orientation and the descriptor for every interest point |
---|
| 180 | for (unsigned n=0; n<ipts.size(); n++){ |
---|
| 181 | //for (Ipoint *k = ipts; k != NULL; k = k->next){ |
---|
| 182 | // set the current interest point |
---|
| 183 | des.setIpoint(&ipts[n]); |
---|
| 184 | // assign reproducible orientation |
---|
| 185 | des.assignOrientation(); |
---|
| 186 | // make the SURF descriptor |
---|
| 187 | des.makeDescriptor(); |
---|
| 188 | } |
---|
| 189 | |
---|
| 190 | // stop measuring the time, we're all done |
---|
| 191 | osmapping::os_GetTime(&tim2); |
---|
| 192 | |
---|
| 193 | // save the interest points in the output file |
---|
| 194 | saveIpoints(fn, ipts, bVerbose, bLaplacian); |
---|
| 195 | |
---|
| 196 | // print some nice information on the command prompt |
---|
| 197 | if( bVerbose ) |
---|
| 198 | cout << "Detection time: " << osmapping::os_TimeDiff(&tim2, &tim1) << " ms" << endl; |
---|
| 199 | |
---|
| 200 | delete im; |
---|
| 201 | |
---|
| 202 | return 0; |
---|
| 203 | } |
---|
| 204 | |
---|
| 205 | // Save the interest points to a regular ASCII file |
---|
| 206 | void saveIpoints(string sFileName, const vector< Ipoint >& ipts, bool bVerbose, bool bLaplacian) |
---|
| 207 | { |
---|
| 208 | ofstream ipfile(sFileName.c_str()); |
---|
| 209 | if( !ipfile ) { |
---|
| 210 | cerr << "ERROR in loadIpoints(): " |
---|
| 211 | << "Couldn't open file '" << sFileName.c_str() << "'!" << endl; //STS |
---|
| 212 | return; |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | double sc; |
---|
| 216 | unsigned count = ipts.size(); |
---|
| 217 | |
---|
| 218 | // Write the file header |
---|
| 219 | if (bLaplacian) |
---|
| 220 | ipfile << VLength + 1 << endl << count << endl; |
---|
| 221 | else |
---|
| 222 | ipfile << VLength << endl << count << endl; |
---|
| 223 | |
---|
| 224 | // In order to just save the interest points without descriptor, comment |
---|
| 225 | // the above and uncomment the following command. |
---|
| 226 | // ipfile << 1.0 << endl << count << endl; |
---|
| 227 | // Save interest point with descriptor in the format of Krystian Mikolajczyk |
---|
| 228 | // for reasons of comparison with other descriptors. As our interest points |
---|
| 229 | // are circular in any case, we use the second component of the ellipse to |
---|
| 230 | // provide some information about the strength of the interest point. This is |
---|
| 231 | // important for 3D reconstruction as only the strongest interest points are |
---|
| 232 | // considered. Replace the strength with 0.0 in order to perform Krystian's |
---|
| 233 | // comparisons. |
---|
| 234 | for (unsigned n=0; n<ipts.size(); n++){ |
---|
| 235 | // circular regions with diameter 5 x scale |
---|
| 236 | sc = 2.5 * ipts[n].scale; sc*=sc; |
---|
| 237 | ipfile << ipts[n].x /* x-location of the interest point */ |
---|
| 238 | << " " << ipts[n].y /* y-location of the interest point */ |
---|
| 239 | << " " << 1.0/sc /* 1/r^2 */ |
---|
| 240 | << " " << 0.0 //(*ipts)[n]->strength /* 0.0 */ |
---|
| 241 | << " " << 1.0/sc; /* 1/r^2 */ |
---|
| 242 | |
---|
| 243 | if (bLaplacian) |
---|
| 244 | ipfile << " " << ipts[n].laplace; |
---|
| 245 | |
---|
| 246 | // Here comes the descriptor |
---|
| 247 | for (int i = 0; i < VLength; i++) { |
---|
| 248 | ipfile << " " << ipts[n].ivec[i]; |
---|
| 249 | } |
---|
| 250 | ipfile << endl; |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | // Write message to terminal. |
---|
| 254 | if( bVerbose ) |
---|
| 255 | cout << count << " interest points found" << endl; |
---|
| 256 | } |
---|
| 257 | |
---|
| 258 | |
---|
| 259 | // Load the interest points from a regular ASCII file |
---|
| 260 | void loadIpoints(string sFileName, vector< Ipoint >& ipts, bool bVerbose) |
---|
| 261 | { |
---|
| 262 | ifstream ipfile(sFileName.c_str()); |
---|
| 263 | if( !ipfile ) { |
---|
| 264 | cerr << "ERROR in loadIpoints(): " |
---|
| 265 | << "Couldn't open file '" << sFileName.c_str() << "'!" << endl; //STS |
---|
| 266 | return; |
---|
| 267 | } |
---|
| 268 | |
---|
| 269 | // Load the file header |
---|
| 270 | float dummy; |
---|
| 271 | unsigned count; |
---|
| 272 | ipfile >> dummy >> count; |
---|
| 273 | |
---|
| 274 | // create a new interest point vector |
---|
| 275 | ipts.clear(); |
---|
| 276 | ipts.resize(count); |
---|
| 277 | |
---|
| 278 | // Load the interest points in Mikolajczyk's format |
---|
| 279 | for (unsigned n=0; n<count; n++){ |
---|
| 280 | // circular regions with diameter 5 x scale |
---|
| 281 | float x, y, a, b, c; |
---|
| 282 | ipfile >> x >> y >> a >> b >> c; |
---|
| 283 | |
---|
| 284 | float det = sqrt((a-c)*(a-c) + 4.0*b*b); |
---|
| 285 | float e1 = 0.5*(a+c + det); |
---|
| 286 | float e2 = 0.5*(a+c - det); |
---|
| 287 | float l1 = (1.0/sqrt(e1)); |
---|
| 288 | float l2 = (1.0/sqrt(e2)); |
---|
| 289 | float sc = sqrt( l1*l2 ); |
---|
| 290 | |
---|
| 291 | ipts[n].x = x; |
---|
| 292 | ipts[n].y = y; |
---|
| 293 | ipts[n].scale = sc/2.5; |
---|
| 294 | } |
---|
| 295 | |
---|
| 296 | // close the interest point file again |
---|
| 297 | ipfile.close(); |
---|
| 298 | |
---|
| 299 | // Write message to terminal. |
---|
| 300 | if( bVerbose ) |
---|
| 301 | cout << "read in " << count << " interest points." << endl; |
---|
| 302 | } |
---|