1 ////////////////////////////////////////////////////////////////////////////
2 // File: SiftPyramid.cpp
3 // Author: Changchang Wu
4 // Description : Implementation of the SiftPyramid class.
5 //
6 //
7 //
8 // Copyright (c) 2007 University of North Carolina at Chapel Hill
9 // All Rights Reserved
10 //
11 // Permission to use, copy, modify and distribute this software and its
12 // documentation for educational, research and non-profit purposes, without
13 // fee, and without a written agreement is hereby granted, provided that the
14 // above copyright notice and the following paragraph appear in all copies.
15 //
16 // The University of North Carolina at Chapel Hill make no representations
17 // about the suitability of this software for any purpose. It is provided
18 // 'as is' without express or implied warranty.
19 //
20 // Please send BUG REPORTS to ccwu@cs.unc.edu
21 //
22 ////////////////////////////////////////////////////////////////////////////
23
24
25 #include "GL/glew.h"
26 #include <string.h>
27 #include <iostream>
28 #include <iomanip>
29 #include <vector>
30 #include <algorithm>
31 #include <fstream>
32 #include <math.h>
33 using namespace std;
34
35 #include "GlobalUtil.h"
36 #include "SiftPyramid.h"
37 #include "SiftGPU.h"
38
39
40 #ifdef DEBUG_SIFTGPU
41 #include "IL/il.h"
42 #include "direct.h"
43 #include "io.h"
44 #include <sys/stat.h>
45 #endif
46
47
48
RunSIFT(GLTexInput * input)49 void SiftPyramid::RunSIFT(GLTexInput*input)
50 {
51 CleanupBeforeSIFT();
52
53 if(_existing_keypoints & SIFT_SKIP_FILTERING)
54 {
55
56 }else
57 {
58 GlobalUtil::StartTimer("Build Pyramid");
59 BuildPyramid(input);
60 GlobalUtil::StopTimer();
61 _timing[0] = GetElapsedTime();
62 }
63
64
65 if(_existing_keypoints)
66 {
67 //existing keypoint list should at least have the locations and scale
68 GlobalUtil::StartTimer("Upload Feature List");
69 if(!(_existing_keypoints & SIFT_SKIP_FILTERING)) ComputeGradient();
70 GenerateFeatureListTex();
71 GlobalUtil::StopTimer();
72 _timing[2] = GetElapsedTime();
73 }else
74 {
75
76 GlobalUtil::StartTimer("Detect Keypoints");
77 DetectKeypointsEX();
78 GlobalUtil::StopTimer();
79 _timing[1] = GetElapsedTime();
80
81 if(GlobalUtil::_ListGenGPU ==1)
82 {
83 GlobalUtil::StartTimer("Get Feature List");
84 GenerateFeatureList();
85 GlobalUtil::StopTimer();
86
87 }else
88 {
89 GlobalUtil::StartTimer("Transfer Feature List");
90 GenerateFeatureListCPU();
91 GlobalUtil::StopTimer();
92 }
93 LimitFeatureCount(0);
94 _timing[2] = GetElapsedTime();
95 }
96
97
98
99 if(_existing_keypoints& SIFT_SKIP_ORIENTATION)
100 {
101 //use exisitng feature orientation or
102 }else if(GlobalUtil::_MaxOrientation>0)
103 {
104 //some extra tricks are done to handle existing keypoint list
105 GlobalUtil::StartTimer("Feature Orientations");
106 GetFeatureOrientations();
107 GlobalUtil::StopTimer();
108 _timing[3] = GetElapsedTime();
109
110 //for existing keypoint list, only the strongest orientation is kept.
111 if(GlobalUtil::_MaxOrientation >1 && !_existing_keypoints && !GlobalUtil::_FixedOrientation)
112 {
113 GlobalUtil::StartTimer("MultiO Feature List");
114 ReshapeFeatureListCPU();
115 LimitFeatureCount(1);
116 GlobalUtil::StopTimer();
117 _timing[4] = GetElapsedTime();
118 }
119 }else
120 {
121 GlobalUtil::StartTimer("Feature Orientations");
122 GetSimplifiedOrientation();
123 GlobalUtil::StopTimer();
124 _timing[3] = GetElapsedTime();
125 }
126
127 PrepareBuffer();
128
129 if(_existing_keypoints & SIFT_SKIP_ORIENTATION)
130 {
131 //no need to read back feature if all fields of keypoints are already specified
132 }else
133 {
134 GlobalUtil::StartTimer("Download Keypoints");
135 #ifdef NO_DUPLICATE_DOWNLOAD
136 if(GlobalUtil::_MaxOrientation < 2 || GlobalUtil::_FixedOrientation)
137 #endif
138 DownloadKeypoints();
139 GlobalUtil::StopTimer();
140 _timing[5] = GetElapsedTime();
141 }
142
143
144
145 if(GlobalUtil::_DescriptorPPT)
146 {
147 //desciprotrs are downloaded in descriptor computation of each level
148 GlobalUtil::StartTimer("Get Descriptor");
149 GetFeatureDescriptors();
150 GlobalUtil::StopTimer();
151 _timing[6] = GetElapsedTime();
152 }
153
154 //reset the existing keypoints
155 _existing_keypoints = 0;
156 _keypoint_index.resize(0);
157
158 if(GlobalUtil::_UseSiftGPUEX)
159 {
160 GlobalUtil::StartTimer("Gen. Display VBO");
161 GenerateFeatureDisplayVBO();
162 GlobalUtil::StopTimer();
163 _timing[7] = GlobalUtil::GetElapsedTime();
164 }
165 //clean up
166 CleanUpAfterSIFT();
167 }
168
169
LimitFeatureCount(int have_keylist)170 void SiftPyramid::LimitFeatureCount(int have_keylist)
171 {
172
173 if(GlobalUtil::_FeatureCountThreshold <= 0 || _existing_keypoints) return;
174 ///////////////////////////////////////////////////////////////
175 //skip the lowest levels to reduce number of features.
176
177 if(GlobalUtil::_TruncateMethod == 2)
178 {
179 int i = 0, new_feature_num = 0, level_num = param._dog_level_num * _octave_num;
180 for(; new_feature_num < _FeatureCountThreshold && i < level_num; ++i) new_feature_num += _levelFeatureNum[i];
181 for(; i < level_num; ++i) _levelFeatureNum[i] = 0;
182
183 if(new_feature_num < _featureNum)
184 {
185 _featureNum = new_feature_num;
186 if(GlobalUtil::_verbose )
187 {
188 std::cout<<"#Features Reduced:\t"<<_featureNum<<endl;
189 }
190 }
191 }else
192 {
193 int i = 0, num_to_erase = 0;
194 while(_featureNum - _levelFeatureNum[i] > _FeatureCountThreshold)
195 {
196 num_to_erase += _levelFeatureNum[i];
197 _featureNum -= _levelFeatureNum[i];
198 _levelFeatureNum[i++] = 0;
199 }
200 if(num_to_erase > 0 && have_keylist)
201 {
202 _keypoint_buffer.erase(_keypoint_buffer.begin(), _keypoint_buffer.begin() + num_to_erase * 4);
203 }
204 if(GlobalUtil::_verbose && num_to_erase > 0)
205 {
206 std::cout<<"#Features Reduced:\t"<<_featureNum<<endl;
207 }
208 }
209
210
211 }
212
PrepareBuffer()213 void SiftPyramid::PrepareBuffer()
214 {
215 //when there is no existing keypoint list, the feature list need to be downloaded
216 //when an existing keypoint list does not have orientaiton, we need to download them again.
217 if(!(_existing_keypoints & SIFT_SKIP_ORIENTATION))
218 {
219 //_keypoint_buffer.resize(4 * (_featureNum +align));
220 _keypoint_buffer.resize(4 * (_featureNum + GlobalUtil::_texMaxDim)); //11/19/2008
221 }
222 if(GlobalUtil::_DescriptorPPT)
223 {
224 //_descriptor_buffer.resize(128*(_featureNum + align));
225 _descriptor_buffer.resize(128 * _featureNum + 16 * GlobalUtil::_texMaxDim);//11/19/2008
226 }
227
228 }
229
GetRequiredOctaveNum(int inputsz)230 int SiftPyramid:: GetRequiredOctaveNum(int inputsz)
231 {
232 //[2 ^ i, 2 ^ (i + 1)) -> i - 3...
233 //768 in [2^9, 2^10) -> 6 -> smallest will be 768 / 32 = 24
234 int num = (int) floor (log ( inputsz * 2.0 / GlobalUtil::_texMinDim )/log(2.0));
235 return num <= 0 ? 1 : num;
236 }
237
CopyFeatureVector(float * keys,float * descriptors)238 void SiftPyramid::CopyFeatureVector(float*keys, float *descriptors)
239 {
240 if(keys) memcpy(keys, &_keypoint_buffer[0], 4*_featureNum*sizeof(float));
241 if(descriptors) memcpy(descriptors, &_descriptor_buffer[0], 128*_featureNum*sizeof(float));
242 }
243
SetKeypointList(int num,const float * keys,int run_on_current,int skip_orientation)244 void SiftPyramid:: SetKeypointList(int num, const float * keys, int run_on_current, int skip_orientation)
245 {
246 //for each input keypoint
247 //sort the key point list by size, and assign them to corresponding levels
248 if(num <=0) return;
249 _featureNum = num;
250 ///copy the keypoints
251 _keypoint_buffer.resize(num * 4);
252 memcpy(&_keypoint_buffer[0], keys, 4 * num * sizeof(float));
253 //location and scale can be skipped
254 _existing_keypoints = SIFT_SKIP_DETECTION;
255 //filtering is skipped if it is running on the same image
256 if(run_on_current) _existing_keypoints |= SIFT_SKIP_FILTERING;
257 //orientation can be skipped if specified
258 if(skip_orientation) _existing_keypoints |= SIFT_SKIP_ORIENTATION;
259 //hacking parameter for using rectangle description mode
260 if(skip_orientation == -1) _existing_keypoints |= SIFT_RECT_DESCRIPTION;
261 }
262
263
SaveSIFT(const char * szFileName)264 void SiftPyramid::SaveSIFT(const char * szFileName)
265 {
266 if (_featureNum <=0) return;
267 float * pk = &_keypoint_buffer[0];
268
269 if(GlobalUtil::_BinarySIFT)
270 {
271 std::ofstream out(szFileName, ios::binary);
272 out.write((char* )(&_featureNum), sizeof(int));
273
274 if(GlobalUtil::_DescriptorPPT)
275 {
276 int dim = 128;
277 out.write((char* )(&dim), sizeof(int));
278 float * pd = &_descriptor_buffer[0] ;
279 for(int i = 0; i < _featureNum; i++, pk+=4, pd +=128)
280 {
281 out.write((char* )(pk +1), sizeof(float));
282 out.write((char* )(pk), sizeof(float));
283 out.write((char* )(pk+2), 2 * sizeof(float));
284 out.write((char* )(pd), 128 * sizeof(float));
285 }
286 }else
287 {
288 int dim = 0;
289 out.write((char* )(&dim), sizeof(int));
290 for(int i = 0; i < _featureNum; i++, pk+=4)
291 {
292 out.write((char* )(pk +1), sizeof(float));
293 out.write((char* )(pk), sizeof(float));
294 out.write((char* )(pk+2), 2 * sizeof(float));
295 }
296 }
297 }else
298 {
299 std::ofstream out(szFileName);
300 out.flags(ios::fixed);
301
302 if(GlobalUtil::_DescriptorPPT)
303 {
304 float * pd = &_descriptor_buffer[0] ;
305 out<<_featureNum<<" 128"<<endl;
306
307 for(int i = 0; i < _featureNum; i++)
308 {
309 //in y, x, scale, orientation order
310 out<<setprecision(2) << pk[1]<<" "<<setprecision(2) << pk[0]<<" "
311 <<setprecision(3) << pk[2]<<" " <<setprecision(3) << pk[3]<< endl;
312
313 ////out << setprecision(12) << pk[1] << " " << pk[0] << " " << pk[2] << " " << pk[3] << endl;
314 pk+=4;
315 for(int k = 0; k < 128; k ++, pd++)
316 {
317 if(GlobalUtil::_NormalizedSIFT)
318 out<< ((unsigned int)floor(0.5+512.0f*(*pd)))<<" ";
319 else
320 out << setprecision(8) << pd[0] << " ";
321
322 if ( (k+1)%20 == 0 ) out<<endl; //suggested by Martin Schneider
323
324 }
325 out<<endl;
326
327 }
328
329 }else
330 {
331 out<<_featureNum<<" 0"<<endl;
332 for(int i = 0; i < _featureNum; i++, pk+=4)
333 {
334 out<<pk[1]<<" "<<pk[0]<<" "<<pk[2]<<" " << pk[3]<<endl;
335 }
336 }
337 }
338 }
339
340 #ifdef DEBUG_SIFTGPU
BeginDEBUG(const char * imagepath)341 void SiftPyramid::BeginDEBUG(const char *imagepath)
342 {
343 if(imagepath && imagepath[0])
344 {
345 strcpy(_debug_path, imagepath);
346 strcat(_debug_path, ".debug");
347 }else
348 {
349 strcpy(_debug_path, ".debug");
350 }
351
352 mkdir(_debug_path);
353 chmod(_debug_path, _S_IREAD | _S_IWRITE);
354 }
355
StopDEBUG()356 void SiftPyramid::StopDEBUG()
357 {
358 _debug_path[0] = 0;
359 }
360
361
WriteTextureForDEBUG(GLTexImage * tex,const char * namet,...)362 void SiftPyramid::WriteTextureForDEBUG(GLTexImage * tex, const char *namet, ...)
363 {
364 char name[_MAX_PATH];
365 char * p = name, * ps = _debug_path;
366 while(*ps) *p++ = *ps ++;
367 *p++ = '/';
368 va_list marker;
369 va_start(marker, namet);
370 vsprintf(p, namet, marker);
371 va_end(marker);
372 unsigned int imID;
373 int width = tex->GetImgWidth();
374 int height = tex->GetImgHeight();
375 float* buffer1 = new float[ width * height * 4];
376 float* buffer2 = new float[ width * height * 4];
377
378 //read data back
379 glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
380 tex->AttachToFBO(0);
381 tex->FitTexViewPort();
382 glReadPixels(0, 0, width, height, GL_RGBA , GL_FLOAT, buffer1);
383
384 //Tiffs saved with IL are flipped
385 for(int i = 0; i < height; i++)
386 {
387 memcpy(buffer2 + i * width * 4,
388 buffer1 + (height - i - 1) * width * 4,
389 width * 4 * sizeof(float));
390 }
391
392 //save data as floating point tiff file
393 ilGenImages(1, &imID);
394 ilBindImage(imID);
395 ilEnable(IL_FILE_OVERWRITE);
396 ilTexImage(width, height, 1, 4, IL_RGBA, IL_FLOAT, buffer2);
397 ilSave(IL_TIF, name);
398 ilDeleteImages(1, &imID);
399
400 delete buffer1;
401 delete buffer2;
402 glReadBuffer(GL_NONE);
403 }
404
405
406 #endif
407