PyramidCU.cpp
Go to the documentation of this file.
00001 
00002 //      File:           PyramidCU.cpp
00003 //      Author:         Changchang Wu
00004 //      Description : implementation of the PyramidCU class.
00005 //                              CUDA-based implementation of SiftPyramid
00006 //
00007 //      Copyright (c) 2007 University of North Carolina at Chapel Hill
00008 //      All Rights Reserved
00009 //
00010 //      Permission to use, copy, modify and distribute this software and its
00011 //      documentation for educational, research and non-profit purposes, without
00012 //      fee, and without a written agreement is hereby granted, provided that the
00013 //      above copyright notice and the following paragraph appear in all copies.
00014 //      
00015 //      The University of North Carolina at Chapel Hill make no representations
00016 //      about the suitability of this software for any purpose. It is provided
00017 //      'as is' without express or implied warranty. 
00018 //
00019 //      Please send BUG REPORTS to ccwu@cs.unc.edu
00020 //
00022 
00023 #if defined(CUDA_SIFTGPU_ENABLED)
00024 
00025 
00026 #include "GL/glew.h"
00027 #include <iostream>
00028 #include <vector>
00029 #include <algorithm>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 #include <math.h>
00033 using namespace std;
00034 
00035 #include "GlobalUtil.h"
00036 #include "GLTexImage.h"
00037 #include "CuTexImage.h" 
00038 #include "SiftGPU.h"
00039 #include "SiftPyramid.h"
00040 #include "ProgramCU.h"
00041 #include "PyramidCU.h"
00042 
00043 
00044 //#include "imdebug/imdebuggl.h"
00045 //#pragma comment (lib, "../lib/imdebug.lib")
00046 
00047 
00048 
00049 #define USE_TIMING()            double t, t0, tt;
00050 #define OCTAVE_START()          if(GlobalUtil::_timingO){       t = t0 = CLOCK();       cout<<"#"<<i+_down_sample_factor<<"\t"; }
00051 #define LEVEL_FINISH()          if(GlobalUtil::_timingL){       ProgramCU::FinishCUDA();        tt = CLOCK();cout<<(tt-t)<<"\t";        t = CLOCK();}
00052 #define OCTAVE_FINISH()         if(GlobalUtil::_timingO)cout<<"|\t"<<(CLOCK()-t0)<<endl;
00053 
00054 
00055 PyramidCU::PyramidCU(SiftParam& sp) : SiftPyramid(sp)
00056 {
00057         _allPyramid = NULL;
00058         _histoPyramidTex = NULL;
00059         _featureTex = NULL;
00060         _descriptorTex = NULL;
00061         _orientationTex = NULL;
00062         _bufferPBO = 0;
00063     _bufferTEX = NULL;
00064         _inputTex = new CuTexImage();
00065 
00067     InitializeContext();
00068 }
00069 
00070 PyramidCU::~PyramidCU()
00071 {
00072         DestroyPerLevelData();
00073         DestroySharedData();
00074         DestroyPyramidData();
00075         if(_inputTex) delete _inputTex;
00076     if(_bufferPBO) glDeleteBuffers(1, &_bufferPBO);
00077     if(_bufferTEX) delete _bufferTEX;
00078 }
00079 
00080 void PyramidCU::InitializeContext()
00081 {
00082     GlobalUtil::InitGLParam(1);
00083     GlobalUtil::_GoodOpenGL = max(GlobalUtil::_GoodOpenGL, 1); 
00084 }
00085 
00086 void PyramidCU::InitPyramid(int w, int h, int ds)
00087 {
00088         int wp, hp, toobig = 0;
00089         if(ds == 0)
00090         {
00091                 //
00092                 TruncateWidth(w);
00094                 _down_sample_factor = 0;
00095                 if(GlobalUtil::_octave_min_default>=0)
00096                 {
00097                         wp = w >> _octave_min_default;
00098                         hp = h >> _octave_min_default;
00099                 }else
00100                 {
00101                         //can't upsample by more than 8
00102                         _octave_min_default = max(-3, _octave_min_default);
00103                         //
00104                         wp = w << (-_octave_min_default);
00105                         hp = h << (-_octave_min_default);
00106                 }
00107                 _octave_min = _octave_min_default;
00108         }else
00109         {
00110                 //must use 0 as _octave_min; 
00111                 _octave_min = 0;
00112                 _down_sample_factor = ds;
00113                 w >>= ds;
00114                 h >>= ds;
00116 
00117                 TruncateWidth(w);
00118 
00119                 wp = w;
00120                 hp = h; 
00121 
00122         }
00123 
00124         while(wp > GlobalUtil::_texMaxDim  || hp > GlobalUtil::_texMaxDim )
00125         {
00126                 _octave_min ++;
00127                 wp >>= 1;
00128                 hp >>= 1;
00129                 toobig = 1;
00130         }
00131         if(toobig && GlobalUtil::_verbose && _octave_min > 0)
00132         {
00133                 std::cout<< "**************************************************************\n"
00134                                         "Image larger than allowed dimension, data will be downsampled!\n"
00135                                         "use -maxd to change the settings\n"
00136                                         "***************************************************************\n";
00137         }
00138         //ResizePyramid(wp, hp);
00139         if( wp == _pyramid_width && hp == _pyramid_height && _allocated )
00140         {
00141                 FitPyramid(wp, hp);
00142         }else if(GlobalUtil::_ForceTightPyramid || _allocated ==0)
00143         {
00144                 ResizePyramid(wp, hp);
00145         }
00146         else if( wp > _pyramid_width || hp > _pyramid_height )
00147         {
00148                 ResizePyramid(max(wp, _pyramid_width), max(hp, _pyramid_height));
00149                 if(wp < _pyramid_width || hp < _pyramid_height)  FitPyramid(wp, hp);
00150         }
00151         else
00152         {
00153                 //try use the pyramid allocated for large image on small input images
00154                 FitPyramid(wp, hp);
00155         }
00156 }
00157 
00158 void PyramidCU::ResizePyramid(int w, int h)
00159 {
00160         //
00161         unsigned int totalkb = 0;
00162         int _octave_num_new, input_sz, i, j;
00163         //
00164 
00165         if(_pyramid_width == w && _pyramid_height == h && _allocated) return;
00166 
00167         if(w > GlobalUtil::_texMaxDim || h > GlobalUtil::_texMaxDim) return ;
00168 
00169         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<w<<"x"<<h<<endl;
00170         //first octave does not change
00171         _pyramid_octave_first = 0;
00172 
00173         
00174         //compute # of octaves
00175 
00176         input_sz = min(w,h) ;
00177         _pyramid_width =  w;
00178         _pyramid_height =  h;
00179 
00180 
00181 
00182         //reset to preset parameters
00183 
00184         _octave_num_new  = GlobalUtil::_octave_num_default;
00185 
00186         if(_octave_num_new < 1) 
00187         {
00188                 _octave_num_new = (int) floor (log ( double(input_sz))/log(2.0)) -3 ;
00189                 if(_octave_num_new<1 ) _octave_num_new = 1;
00190         }
00191 
00192         if(_pyramid_octave_num != _octave_num_new)
00193         {
00194                 //destroy the original pyramid if the # of octave changes
00195                 if(_octave_num >0) 
00196                 {
00197                         DestroyPerLevelData();
00198                         DestroyPyramidData();
00199                 }
00200                 _pyramid_octave_num = _octave_num_new;
00201         }
00202 
00203         _octave_num = _pyramid_octave_num;
00204 
00205         int noct = _octave_num;
00206         int nlev = param._level_num;
00207 
00208         //      //initialize the pyramid
00209         if(_allPyramid==NULL)   _allPyramid = new CuTexImage[ noct* nlev * DATA_NUM];
00210 
00211         CuTexImage * gus =  GetBaseLevel(_octave_min, DATA_GAUSSIAN);
00212         CuTexImage * dog =  GetBaseLevel(_octave_min, DATA_DOG);
00213         CuTexImage * got =  GetBaseLevel(_octave_min, DATA_GRAD);
00214         CuTexImage * key =  GetBaseLevel(_octave_min, DATA_KEYPOINT);
00215 
00217 
00218         for(i = 0; i< noct; i++)
00219         {
00220                 int wa = ((w + 3) / 4) * 4;
00221 
00222                 totalkb += ((nlev *8 -19)* (wa * h) * 4 / 1024);
00223                 for( j = 0; j< nlev; j++, gus++, dog++, got++, key++)
00224                 {
00225                         gus->InitTexture(wa, h); //nlev
00226                         if(j==0)continue;
00227                         dog->InitTexture(wa, h);  //nlev -1
00228                         if(     j >= 1 && j < 1 + param._dog_level_num)
00229                         {
00230                                 got->InitTexture(wa, h, 2); //2 * nlev - 6
00231                                 got->InitTexture2D();
00232                         }
00233                         if(j > 1 && j < nlev -1)        key->InitTexture(wa, h, 4); // nlev -3 ; 4 * nlev - 12
00234                 }
00235                 w>>=1;
00236                 h>>=1;
00237         }
00238 
00239         totalkb += ResizeFeatureStorage();
00240 
00241         if(ProgramCU::CheckErrorCUDA("ResizePyramid")) SetFailStatus(); 
00242 
00243         _allocated = 1;
00244 
00245         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<(totalkb/1024)<<"MB\n";
00246 
00247 }
00248 
00249 void PyramidCU::FitPyramid(int w, int h)
00250 {
00251         _pyramid_octave_first = 0;
00252         //
00253         _octave_num  = GlobalUtil::_octave_num_default;
00254 
00255         int _octave_num_max = max(1, (int) floor (log ( double(min(w, h)))/log(2.0))  -3 );
00256 
00257         if(_octave_num < 1 || _octave_num > _octave_num_max) 
00258         {
00259                 _octave_num = _octave_num_max;
00260         }
00261 
00262 
00263         int pw = _pyramid_width>>1, ph = _pyramid_height>>1;
00264         while(_pyramid_octave_first + _octave_num < _pyramid_octave_num &&  
00265                 pw >= w && ph >= h)
00266         {
00267                 _pyramid_octave_first++;
00268                 pw >>= 1;
00269                 ph >>= 1;
00270         }
00271 
00273         int nlev = param._level_num;
00274         CuTexImage * gus =  GetBaseLevel(_octave_min, DATA_GAUSSIAN);
00275         CuTexImage * dog =  GetBaseLevel(_octave_min, DATA_DOG);
00276         CuTexImage * got =  GetBaseLevel(_octave_min, DATA_GRAD);
00277         CuTexImage * key =  GetBaseLevel(_octave_min, DATA_KEYPOINT);
00278         for(int i = 0; i< _octave_num; i++)
00279         {
00280                 int wa = ((w + 3) / 4) * 4;
00281 
00282                 for(int j = 0; j< nlev; j++, gus++, dog++, got++, key++)
00283                 {
00284                         gus->InitTexture(wa, h); //nlev
00285                         if(j==0)continue;
00286                         dog->InitTexture(wa, h);  //nlev -1
00287                         if(     j >= 1 && j < 1 + param._dog_level_num)
00288                         {
00289                                 got->InitTexture(wa, h, 2); //2 * nlev - 6
00290                                 got->InitTexture2D();
00291                         }
00292                         if(j > 1 && j < nlev -1)        key->InitTexture(wa, h, 4); // nlev -3 ; 4 * nlev - 12
00293                 }
00294                 w>>=1;
00295                 h>>=1;
00296         }
00297 }
00298 
00299 int PyramidCU::CheckCudaDevice(int device)
00300 {
00301     return ProgramCU::CheckCudaDevice(device);
00302 }
00303 
00304 void PyramidCU::SetLevelFeatureNum(int idx, int fcount)
00305 {
00306         _featureTex[idx].InitTexture(fcount, 1, 4);
00307         _levelFeatureNum[idx] = fcount;
00308 }
00309 
00310 int PyramidCU::ResizeFeatureStorage()
00311 {
00312         int totalkb = 0;
00313         if(_levelFeatureNum==NULL)      _levelFeatureNum = new int[_octave_num * param._dog_level_num];
00314         std::fill(_levelFeatureNum, _levelFeatureNum+_octave_num * param._dog_level_num, 0); 
00315 
00316         int wmax = GetBaseLevel(_octave_min)->GetImgWidth();
00317         int hmax = GetBaseLevel(_octave_min)->GetImgHeight();
00318         int whmax = max(wmax, hmax);
00319         int w,  i;
00320 
00321         //
00322         int num = (int)ceil(log(double(whmax))/log(4.0));
00323 
00324         if( _hpLevelNum != num)
00325         {
00326                 _hpLevelNum = num;
00327                 if(_histoPyramidTex ) delete [] _histoPyramidTex;
00328                 _histoPyramidTex = new CuTexImage[_hpLevelNum];
00329         }
00330 
00331         for(i = 0, w = 1; i < _hpLevelNum; i++)
00332         {
00333                 _histoPyramidTex[i].InitTexture(w, whmax, 4);
00334                 w<<=2;
00335         }
00336 
00337         // (4 ^ (_hpLevelNum) -1 / 3) pixels
00338         totalkb += (((1 << (2 * _hpLevelNum)) -1) / 3 * 16 / 1024);
00339 
00340         //initialize the feature texture
00341         int idx = 0, n = _octave_num * param._dog_level_num;
00342         if(_featureTex==NULL)   _featureTex = new CuTexImage[n];
00343         if(GlobalUtil::_MaxOrientation >1 && GlobalUtil::_OrientationPack2==0 && _orientationTex== NULL)
00344                 _orientationTex = new CuTexImage[n];
00345 
00346 
00347         for(i = 0; i < _octave_num; i++)
00348         {
00349                 CuTexImage * tex = GetBaseLevel(i+_octave_min);
00350                 int fmax = int(tex->GetImgWidth() * tex->GetImgHeight()*GlobalUtil::_MaxFeaturePercent);
00351                 //
00352                 if(fmax > GlobalUtil::_MaxLevelFeatureNum) fmax = GlobalUtil::_MaxLevelFeatureNum;
00353                 else if(fmax < 32) fmax = 32;   //give it at least a space of 32 feature
00354 
00355                 for(int j = 0; j < param._dog_level_num; j++, idx++)
00356                 {
00357                         _featureTex[idx].InitTexture(fmax, 1, 4);
00358                         totalkb += fmax * 16 /1024;
00359                         //
00360                         if(GlobalUtil::_MaxOrientation>1 && GlobalUtil::_OrientationPack2 == 0)
00361                         {
00362                                 _orientationTex[idx].InitTexture(fmax, 1, 4);
00363                                 totalkb += fmax * 16 /1024;
00364                         }
00365                 }
00366         }
00367 
00368 
00369         //this just need be initialized once
00370         if(_descriptorTex==NULL)
00371         {
00372                 //initialize feature texture pyramid
00373                 int fmax = _featureTex->GetImgWidth();
00374                 _descriptorTex = new CuTexImage;
00375                 totalkb += ( fmax /2);
00376                 _descriptorTex->InitTexture(fmax *128, 1, 1);
00377         }else
00378         {
00379                 totalkb +=  _descriptorTex->GetDataSize()/1024;
00380         }
00381         return totalkb;
00382 }
00383 
00384 void PyramidCU::GetFeatureDescriptors() 
00385 {
00386         //descriptors...
00387         float* pd =  &_descriptor_buffer[0];
00388         vector<float> descriptor_buffer2;
00389 
00390         //use another buffer if we need to re-order the descriptors
00391         if(_keypoint_index.size() > 0)
00392         {
00393                 descriptor_buffer2.resize(_descriptor_buffer.size());
00394                 pd = &descriptor_buffer2[0];
00395         }
00396 
00397         CuTexImage * got, * ftex= _featureTex;
00398         for(int i = 0, idx = 0; i < _octave_num; i++)
00399         {
00400                 got = GetBaseLevel(i + _octave_min, DATA_GRAD) + 1;
00401                 for(int j = 0; j < param._dog_level_num; j++, ftex++, idx++, got++)
00402                 {
00403                         if(_levelFeatureNum[idx]==0) continue;
00404             ProgramCU::ComputeDescriptor(ftex, got, _descriptorTex, IsUsingRectDescription());//process
00405                         _descriptorTex->CopyToHost(pd); //readback descriptor
00406                         pd += 128*_levelFeatureNum[idx];
00407                 }
00408         }
00409 
00410         if(GlobalUtil::_timingS) ProgramCU::FinishCUDA();
00411 
00412         if(_keypoint_index.size() > 0)
00413         {
00414             //put the descriptor back to the original order for keypoint list.
00415                 for(int i = 0; i < _featureNum; ++i)
00416                 {
00417                         int index = _keypoint_index[i];
00418                         memcpy(&_descriptor_buffer[index*128], &descriptor_buffer2[i*128], 128 * sizeof(float));
00419                 }
00420         }
00421 
00422         if(ProgramCU::CheckErrorCUDA("PyramidCU::GetFeatureDescriptors")) SetFailStatus(); 
00423 }
00424 
00425 void PyramidCU::GenerateFeatureListTex() 
00426 {
00427 
00428         vector<float> list;
00429         int idx = 0;
00430         const double twopi = 2.0*3.14159265358979323846;
00431         float sigma_half_step = powf(2.0f, 0.5f / param._dog_level_num);
00432         float octave_sigma = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
00433         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f; 
00434         if(_down_sample_factor>0) octave_sigma *= float(1<<_down_sample_factor); 
00435 
00436         _keypoint_index.resize(0); // should already be 0
00437         for(int i = 0; i < _octave_num; i++, octave_sigma*= 2.0f)
00438         {
00439                 for(int j = 0; j < param._dog_level_num; j++, idx++)
00440                 {
00441                         list.resize(0);
00442                         float level_sigma = param.GetLevelSigma(j + param._level_min + 1) * octave_sigma;
00443                         float sigma_min = level_sigma / sigma_half_step;
00444                         float sigma_max = level_sigma * sigma_half_step;
00445                         int fcount = 0 ;
00446                         for(int k = 0; k < _featureNum; k++)
00447                         {
00448                                 float * key = &_keypoint_buffer[k*4];
00449                 float sigmak = key[2]; 
00451                 if(IsUsingRectDescription()) sigmak = min(key[2], key[3]) / 12.0f; 
00452 
00453                                 if(   (sigmak >= sigma_min && sigmak < sigma_max)
00454                                         ||(sigmak < sigma_min && i ==0 && j == 0)
00455                                         ||(sigmak > sigma_max && i == _octave_num -1 && j == param._dog_level_num - 1))
00456                                 {
00457                                         //add this keypoint to the list
00458                                         list.push_back((key[0] - offset) / octave_sigma + 0.5f);
00459                                         list.push_back((key[1] - offset) / octave_sigma + 0.5f);
00460                     if(IsUsingRectDescription())
00461                     {
00462                         list.push_back(key[2] / octave_sigma);
00463                         list.push_back(key[3] / octave_sigma);
00464                     }else
00465                     {
00466                                             list.push_back(key[2] / octave_sigma);
00467                                             list.push_back((float)fmod(twopi-key[3], twopi));
00468                     }
00469                                         fcount ++;
00470                                         //save the index of keypoints
00471                                         _keypoint_index.push_back(k);
00472                                 }
00473 
00474                         }
00475 
00476                         _levelFeatureNum[idx] = fcount;
00477                         if(fcount==0)continue;
00478                         CuTexImage * ftex = _featureTex+idx;
00479 
00480                         SetLevelFeatureNum(idx, fcount);
00481                         ftex->CopyFromHost(&list[0]);
00482                 }
00483         }
00484 
00485         if(GlobalUtil::_verbose)
00486         {
00487                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
00488         }
00489 
00490 }
00491 
00492 void PyramidCU::ReshapeFeatureListCPU() 
00493 {
00494         int i, szmax =0, sz;
00495         int n = param._dog_level_num*_octave_num;
00496         for( i = 0; i < n; i++) 
00497         {
00498                 sz = _levelFeatureNum[i];
00499                 if(sz > szmax ) szmax = sz;
00500         }
00501         float * buffer = new float[szmax*16];
00502         float * buffer1 = buffer;
00503         float * buffer2 = buffer + szmax*4;
00504 
00505 
00506 
00507         _featureNum = 0;
00508 
00509 #ifdef NO_DUPLICATE_DOWNLOAD
00510         const double twopi = 2.0*3.14159265358979323846;
00511         _keypoint_buffer.resize(0);
00512         float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
00513         if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
00514         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
00515 #endif
00516 
00517 
00518         for(i = 0; i < n; i++)
00519         {
00520                 if(_levelFeatureNum[i]==0)continue;
00521 
00522                 _featureTex[i].CopyToHost(buffer1);
00523                 
00524                 int fcount =0;
00525                 float * src = buffer1;
00526                 float * des = buffer2;
00527                 const static double factor  = 2.0*3.14159265358979323846/65535.0;
00528                 for(int j = 0; j < _levelFeatureNum[i]; j++, src+=4)
00529                 {
00530                         unsigned short * orientations = (unsigned short*) (&src[3]);
00531                         if(orientations[0] != 65535)
00532                         {
00533                                 des[0] = src[0];
00534                                 des[1] = src[1];
00535                                 des[2] = src[2];
00536                                 des[3] = float( factor* orientations[0]);
00537                                 fcount++;
00538                                 des += 4;
00539                                 if(orientations[1] != 65535 && orientations[1] != orientations[0])
00540                                 {
00541                                         des[0] = src[0];
00542                                         des[1] = src[1];
00543                                         des[2] = src[2];
00544                                         des[3] = float(factor* orientations[1]);        
00545                                         fcount++;
00546                                         des += 4;
00547                                 }
00548                         }
00549                 }
00550                 //texture size
00551                 SetLevelFeatureNum(i, fcount);
00552                 _featureTex[i].CopyFromHost(buffer2);
00553                 
00554                 if(fcount == 0) continue;
00555 
00556 #ifdef NO_DUPLICATE_DOWNLOAD
00557                 float oss = os * (1 << (i / param._dog_level_num));
00558                 _keypoint_buffer.resize((_featureNum + fcount) * 4);
00559                 float* ds = &_keypoint_buffer[_featureNum * 4];
00560                 float* fs = buffer2;
00561                 for(int k = 0;  k < fcount; k++, ds+=4, fs+=4)
00562                 {
00563                         ds[0] = oss*(fs[0]-0.5f) + offset;      //x
00564                         ds[1] = oss*(fs[1]-0.5f) + offset;      //y
00565                         ds[2] = oss*fs[2];  //scale
00566                         ds[3] = (float)fmod(twopi-fs[3], twopi);        //orientation, mirrored
00567                 }
00568 #endif
00569                 _featureNum += fcount;
00570         }
00571         delete[] buffer;
00572         if(GlobalUtil::_verbose)
00573         {
00574                 std::cout<<"#Features MO:\t"<<_featureNum<<endl;
00575         }
00576 }
00577 
00578 void PyramidCU::GenerateFeatureDisplayVBO() 
00579 {
00580         //it is weried that this part is very slow.
00581         //use a big VBO to save all the SIFT box vertices
00582         int nvbo = _octave_num * param._dog_level_num;
00583         if(_featureDisplayVBO==NULL)
00584         {
00585                 //initialize the vbos
00586                 _featureDisplayVBO = new GLuint[nvbo];
00587                 _featurePointVBO = new GLuint[nvbo];
00588                 glGenBuffers(nvbo, _featureDisplayVBO); 
00589                 glGenBuffers(nvbo, _featurePointVBO);
00590         }
00591         for(int i = 0; i < nvbo; i++)
00592         {
00593                 if(_levelFeatureNum[i]<=0)continue;
00594                 CuTexImage * ftex  = _featureTex + i;
00595                 CuTexImage texPBO1( _levelFeatureNum[i]* 10, 1, 4, _featureDisplayVBO[i]);
00596                 CuTexImage texPBO2(_levelFeatureNum[i], 1, 4, _featurePointVBO[i]);
00597                 ProgramCU::DisplayKeyBox(ftex, &texPBO1);
00598                 ProgramCU::DisplayKeyPoint(ftex, &texPBO2);     
00599         }
00600 }
00601 
00602 void PyramidCU::DestroySharedData() 
00603 {
00604         //histogram reduction
00605         if(_histoPyramidTex)
00606         {
00607                 delete[]        _histoPyramidTex;
00608                 _hpLevelNum = 0;
00609                 _histoPyramidTex = NULL;
00610         }
00611         //descriptor storage shared by all levels
00612         if(_descriptorTex)
00613         {
00614                 delete _descriptorTex;
00615                 _descriptorTex = NULL;
00616         }
00617         //cpu reduction buffer.
00618         if(_histo_buffer)
00619         {
00620                 delete[] _histo_buffer;
00621                 _histo_buffer = 0;
00622         }
00623 }
00624 
00625 void PyramidCU::DestroyPerLevelData() 
00626 {
00627         //integers vector to store the feature numbers.
00628         if(_levelFeatureNum)
00629         {
00630                 delete [] _levelFeatureNum;
00631                 _levelFeatureNum = NULL;
00632         }
00633         //texture used to store features
00634         if(     _featureTex)
00635         {
00636                 delete [] _featureTex;
00637                 _featureTex =   NULL;
00638         }
00639         //texture used for multi-orientation 
00640         if(_orientationTex)
00641         {
00642                 delete [] _orientationTex;
00643                 _orientationTex = NULL;
00644         }
00645         int no = _octave_num* param._dog_level_num;
00646 
00647         //two sets of vbos used to display the features
00648         if(_featureDisplayVBO)
00649         {
00650                 glDeleteBuffers(no, _featureDisplayVBO);
00651                 delete [] _featureDisplayVBO;
00652                 _featureDisplayVBO = NULL;
00653         }
00654         if( _featurePointVBO)
00655         {
00656                 glDeleteBuffers(no, _featurePointVBO);
00657                 delete [] _featurePointVBO;
00658                 _featurePointVBO = NULL;
00659         }
00660 }
00661 
00662 void PyramidCU::DestroyPyramidData()
00663 {
00664         if(_allPyramid)
00665         {
00666                 delete [] _allPyramid;
00667                 _allPyramid = NULL;
00668         }
00669 }
00670 
00671 void PyramidCU::DownloadKeypoints() 
00672 {
00673         const double twopi = 2.0*3.14159265358979323846;
00674         int idx = 0;
00675         float * buffer = &_keypoint_buffer[0];
00676         vector<float> keypoint_buffer2;
00677         //use a different keypoint buffer when processing with an exisint features list
00678         //without orientation information. 
00679         if(_keypoint_index.size() > 0)
00680         {
00681                 keypoint_buffer2.resize(_keypoint_buffer.size());
00682                 buffer = &keypoint_buffer2[0];
00683         }
00684         float * p = buffer, *ps;
00685         CuTexImage * ftex = _featureTex;
00687         float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
00688         if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
00689         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
00691         for(int i = 0; i < _octave_num; i++, os *= 2.0f)
00692         {
00693         
00694                 for(int j = 0; j  < param._dog_level_num; j++, idx++, ftex++)
00695                 {
00696 
00697                         if(_levelFeatureNum[idx]>0)
00698                         {       
00699                                 ftex->CopyToHost(ps = p);
00700                                 for(int k = 0;  k < _levelFeatureNum[idx]; k++, ps+=4)
00701                                 {
00702                                         ps[0] = os*(ps[0]-0.5f) + offset;       //x
00703                                         ps[1] = os*(ps[1]-0.5f) + offset;       //y
00704                                         ps[2] = os*ps[2]; 
00705                                         ps[3] = (float)fmod(twopi-ps[3], twopi);        //orientation, mirrored
00706                                 }
00707                                 p+= 4* _levelFeatureNum[idx];
00708                         }
00709                 }
00710         }
00711 
00712         //put the feature into their original order for existing keypoint 
00713         if(_keypoint_index.size() > 0)
00714         {
00715                 for(int i = 0; i < _featureNum; ++i)
00716                 {
00717                         int index = _keypoint_index[i];
00718                         memcpy(&_keypoint_buffer[index*4], &keypoint_buffer2[i*4], 4 * sizeof(float));
00719                 }
00720         }
00721 }
00722 
00723 void PyramidCU::GenerateFeatureListCPU()
00724 {
00725         //no cpu version provided
00726         GenerateFeatureList();
00727 }
00728 
00729 void PyramidCU::GenerateFeatureList(int i, int j, int reduction_count, vector<int>& hbuffer)
00730 {
00731     int fcount = 0, idx = i * param._dog_level_num  + j;
00732         int hist_level_num = _hpLevelNum - _pyramid_octave_first /2; 
00733         int ii, k, len; 
00734 
00735         CuTexImage * htex, * ftex, * tex, *got;
00736         ftex = _featureTex + idx;
00737         htex = _histoPyramidTex + hist_level_num -1;
00738         tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2 + j;
00739         got = GetBaseLevel(_octave_min + i, DATA_GRAD) + 2 + j;
00740 
00741         ProgramCU::InitHistogram(tex, htex);
00742 
00743         for(k = 0; k < reduction_count - 1; k++, htex--)
00744         {
00745                 ProgramCU::ReduceHistogram(htex, htex -1);      
00746         }
00747         
00748         //htex has the row reduction result
00749         len = htex->GetImgHeight() * 4;
00750         hbuffer.resize(len);
00751         ProgramCU::FinishCUDA();
00752         htex->CopyToHost(&hbuffer[0]);
00753         
00755         for(ii = 0; ii < len; ++ii)     {if(!(hbuffer[ii]>= 0)) hbuffer[ii] = 0; }//?
00756         
00757     
00758     for(ii = 0; ii < len; ++ii)         fcount += hbuffer[ii];
00759         SetLevelFeatureNum(idx, fcount);
00760         
00761     //build the feature list
00762         if(fcount > 0)
00763         {
00764                 _featureNum += fcount;
00765                 _keypoint_buffer.resize(fcount * 4);
00766                 //vector<int> ikbuf(fcount*4);
00767                 int* ibuf = (int*) (&_keypoint_buffer[0]);
00768 
00769                 for(ii = 0; ii < len; ++ii)
00770                 {
00771                         int x = ii%4, y = ii / 4;
00772                         for(int jj = 0 ; jj < hbuffer[ii]; ++jj, ibuf+=4)
00773                         {
00774                                 ibuf[0] = x; ibuf[1] = y; ibuf[2] = jj; ibuf[3] = 0;
00775                         }
00776                 }
00777                 _featureTex[idx].CopyFromHost(&_keypoint_buffer[0]);
00778         
00780                 ProgramCU::GenerateList(_featureTex + idx, ++htex);
00781                 for(k = 2; k < reduction_count; k++)
00782                 {
00783                         ProgramCU::GenerateList(_featureTex + idx, ++htex);
00784                 }
00785         }
00786 }
00787 
00788 void PyramidCU::GenerateFeatureList()
00789 {
00790         double t1, t2; 
00791         int ocount = 0, reduction_count;
00792     int reverse = (GlobalUtil::_TruncateMethod == 1);
00793 
00794         vector<int> hbuffer;
00795         _featureNum = 0;
00796 
00797         //for(int i = 0, idx = 0; i < _octave_num; i++)
00798     FOR_EACH_OCTAVE(i, reverse)
00799         {
00800         CuTexImage* tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2;
00801                 reduction_count = FitHistogramPyramid(tex);
00802 
00803                 if(GlobalUtil::_timingO)
00804                 {
00805                         t1 = CLOCK(); 
00806                         ocount = 0;
00807                         std::cout<<"#"<<i+_octave_min + _down_sample_factor<<":\t";
00808                 }
00809                 //for(int j = 0; j < param._dog_level_num; j++, idx++)
00810         FOR_EACH_LEVEL(j, reverse)
00811                 {
00812             if(GlobalUtil::_TruncateMethod && GlobalUtil::_FeatureCountThreshold > 0 && _featureNum > GlobalUtil::_FeatureCountThreshold) continue;
00813 
00814                 GenerateFeatureList(i, j, reduction_count, hbuffer);
00815 
00817                         if(GlobalUtil::_timingO)
00818                         {
00819                 int idx = i * param._dog_level_num + j;
00820                                 ocount += _levelFeatureNum[idx];
00821                                 std::cout<< _levelFeatureNum[idx] <<"\t";
00822                         }
00823                 }
00824                 if(GlobalUtil::_timingO)
00825                 {       
00826                         t2 = CLOCK(); 
00827                         std::cout << "| \t" << int(ocount) << " :\t(" << (t2 - t1) << ")\n";
00828                 }
00829         }
00831         CopyGradientTex();
00833         if(GlobalUtil::_timingS)ProgramCU::FinishCUDA();
00834 
00835         if(GlobalUtil::_verbose)
00836         {
00837                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
00838         }
00839 
00840         if(ProgramCU::CheckErrorCUDA("PyramidCU::GenerateFeatureList")) SetFailStatus();
00841 }
00842 
00843 GLTexImage* PyramidCU::GetLevelTexture(int octave, int level)
00844 {
00845         return GetLevelTexture(octave, level, DATA_GAUSSIAN);
00846 }
00847 
00848 GLTexImage* PyramidCU::ConvertTexCU2GL(CuTexImage* tex, int dataName)
00849 {
00850         
00851         GLenum format = GL_LUMINANCE;
00852         int convert_done = 1;
00853     if(_bufferPBO == 0) glGenBuffers(1, &_bufferPBO);
00854     if(_bufferTEX == NULL) _bufferTEX = new GLTexImage;
00855         switch(dataName)
00856         {
00857         case DATA_GAUSSIAN:
00858                 {
00859                         convert_done = tex->CopyToPBO(_bufferPBO);
00860                         break;
00861                 }
00862         case DATA_DOG:
00863                 {
00864                         CuTexImage texPBO(tex->GetImgWidth(), tex->GetImgHeight(), 1, _bufferPBO);
00865                         if(texPBO._cuData == 0 || tex->_cuData == NULL) convert_done = 0;
00866                         else ProgramCU::DisplayConvertDOG(tex, &texPBO);
00867                         break;
00868                 }
00869         case DATA_GRAD:
00870                 {
00871                         CuTexImage texPBO(tex->GetImgWidth(), tex->GetImgHeight(), 1, _bufferPBO);
00872                         if(texPBO._cuData == 0 || tex->_cuData == NULL) convert_done = 0;
00873                         else ProgramCU::DisplayConvertGRD(tex, &texPBO);
00874                         break;
00875                 }
00876         case DATA_KEYPOINT:
00877                 {
00878                         CuTexImage * dog = tex - param._level_num * _pyramid_octave_num;
00879                         format = GL_RGBA;
00880                         CuTexImage texPBO(tex->GetImgWidth(), tex->GetImgHeight(), 4, _bufferPBO);
00881                         if(texPBO._cuData == 0 || tex->_cuData == NULL) convert_done = 0;
00882                         else ProgramCU::DisplayConvertKEY(tex, dog, &texPBO);
00883                         break;
00884                 }
00885         default:
00886                         convert_done = 0;
00887                         break;
00888         }
00889 
00890         if(convert_done)
00891         {
00892                 _bufferTEX->InitTexture(max(_bufferTEX->GetTexWidth(), tex->GetImgWidth()), max(_bufferTEX->GetTexHeight(), tex->GetImgHeight()));
00893                 _bufferTEX->CopyFromPBO(_bufferPBO, tex->GetImgWidth(), tex->GetImgHeight(), format);
00894         }else
00895         {
00896                 _bufferTEX->SetImageSize(0, 0);
00897         }
00898 
00899         return _bufferTEX;
00900 }
00901 
00902 GLTexImage* PyramidCU::GetLevelTexture(int octave, int level, int dataName) 
00903 {
00904         CuTexImage* tex = GetBaseLevel(octave, dataName) + (level - param._level_min);
00905         //CuTexImage* gus = GetBaseLevel(octave, DATA_GAUSSIAN) + (level - param._level_min); 
00906         return ConvertTexCU2GL(tex, dataName);
00907 }
00908 
00909 void PyramidCU::ConvertInputToCU(GLTexInput* input)
00910 {
00911         int ws = input->GetImgWidth(), hs = input->GetImgHeight();
00912         TruncateWidth(ws);
00913         //copy the input image to pixel buffer object
00914     if(input->_pixel_data)
00915     {
00916         _inputTex->InitTexture(ws, hs, 1);
00917         _inputTex->CopyFromHost(input->_pixel_data); 
00918     }else
00919     {
00920         if(_bufferPBO == 0) glGenBuffers(1, &_bufferPBO);
00921         if(input->_rgb_converted && input->CopyToPBO(_bufferPBO, ws, hs, GL_LUMINANCE))
00922         {
00923                     _inputTex->InitTexture(ws, hs, 1);
00924             _inputTex->CopyFromPBO(ws, hs, _bufferPBO); 
00925         }else if(input->CopyToPBO(_bufferPBO, ws, hs))
00926             {
00927                     CuTexImage texPBO(ws, hs, 4, _bufferPBO);
00928                     _inputTex->InitTexture(ws, hs, 1);
00929                     ProgramCU::ReduceToSingleChannel(_inputTex, &texPBO, !input->_rgb_converted);
00930             }else
00931             {
00932                     std::cerr<< "Unable To Convert Intput\n";
00933             }
00934     }
00935 }
00936 
00937 void PyramidCU::BuildPyramid(GLTexInput * input)
00938 {
00939 
00940         USE_TIMING();
00941 
00942         int i, j;
00943         
00944         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
00945         {
00946 
00947                 float* filter_sigma = param._sigma;
00948                 CuTexImage *tex = GetBaseLevel(i);
00949                 CuTexImage *buf = GetBaseLevel(i, DATA_KEYPOINT) +2;
00950                 j = param._level_min + 1;
00951 
00952                 OCTAVE_START();
00953 
00954                 if( i == _octave_min )
00955                 {       
00956                         ConvertInputToCU(input);
00957 
00958                         if(i == 0)
00959                         {
00960                                 ProgramCU::FilterImage(tex, _inputTex, buf, 
00961                     param.GetInitialSmoothSigma(_octave_min + _down_sample_factor));
00962                         }else
00963                         {
00964                                 if(i < 0)       ProgramCU::SampleImageU(tex, _inputTex, -i);                    
00965                                 else            ProgramCU::SampleImageD(tex, _inputTex, i);
00966                                 ProgramCU::FilterImage(tex, tex, buf, 
00967                     param.GetInitialSmoothSigma(_octave_min + _down_sample_factor));
00968                         }
00969                 }else
00970                 {
00971                         ProgramCU::SampleImageD(tex, GetBaseLevel(i - 1) + param._level_ds - param._level_min); 
00972                         if(param._sigma_skip1 > 0)
00973                         {
00974                                 ProgramCU::FilterImage(tex, tex, buf, param._sigma_skip1);
00975                         }
00976                 }
00977                 LEVEL_FINISH();
00978                 for( ; j <=  param._level_max ; j++, tex++, filter_sigma++)
00979                 {
00980                         // filtering
00981                         ProgramCU::FilterImage(tex + 1, tex, buf, *filter_sigma);
00982                         LEVEL_FINISH();
00983                 }
00984                 OCTAVE_FINISH();
00985         }
00986         if(GlobalUtil::_timingS) ProgramCU::FinishCUDA();
00987 
00988         if(ProgramCU::CheckErrorCUDA("PyramidCU::BuildPyramid")) SetFailStatus();
00989 }
00990 
00991 void PyramidCU::DetectKeypointsEX()
00992 {
00993 
00994 
00995         int i, j;
00996         double t0, t, ts, t1, t2;
00997 
00998         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
00999 
01000         for(i = _octave_min; i < _octave_min + _octave_num; i++)
01001         {
01002                 CuTexImage * gus = GetBaseLevel(i) + 1;
01003                 CuTexImage * dog = GetBaseLevel(i, DATA_DOG) + 1;
01004                 CuTexImage * got = GetBaseLevel(i, DATA_GRAD) + 1;
01005                 //compute the gradient
01006                 for(j = param._level_min +1; j <=  param._level_max ; j++, gus++, dog++, got++)
01007                 {
01008                         //input: gus and gus -1
01009                         //output: gradient, dog, orientation
01010                         ProgramCU::ComputeDOG(gus, dog, got);
01011                 }
01012         }
01013         if(GlobalUtil::_timingS && GlobalUtil::_verbose)
01014         {
01015                 ProgramCU::FinishCUDA();
01016                 t1 = CLOCK();
01017         }
01018 
01019         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
01020         {
01021                 if(GlobalUtil::_timingO)
01022                 {
01023                         t0 = CLOCK();
01024                         std::cout<<"#"<<(i + _down_sample_factor)<<"\t";
01025                 }
01026                 CuTexImage * dog = GetBaseLevel(i, DATA_DOG) + 2;
01027                 CuTexImage * key = GetBaseLevel(i, DATA_KEYPOINT) +2;
01028 
01029 
01030                 for( j = param._level_min +2; j <  param._level_max ; j++, dog++, key++)
01031                 {
01032                         if(GlobalUtil::_timingL)t = CLOCK();
01033                         //input, dog, dog + 1, dog -1
01034                         //output, key
01035                         ProgramCU::ComputeKEY(dog, key, param._dog_threshold, param._edge_threshold);
01036                         if(GlobalUtil::_timingL)
01037                         {
01038                                 std::cout<<(CLOCK()-t)<<"\t";
01039                         }
01040                 }
01041                 if(GlobalUtil::_timingO)
01042                 {
01043                         std::cout<<"|\t"<<(CLOCK()-t0)<<"\n";
01044                 }
01045         }
01046 
01047         if(GlobalUtil::_timingS)
01048         {
01049                 ProgramCU::FinishCUDA();
01050                 if(GlobalUtil::_verbose) 
01051                 {       
01052                         t2 = CLOCK();
01053                         std::cout       <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n"
01054                                                 <<"<Get Keypoints  >\t"<<(t2-t1)<<"\n";
01055                 }                               
01056         }
01057 }
01058 
01059 void PyramidCU::CopyGradientTex()
01060 {
01061         double ts, t1;
01062 
01063         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
01064 
01065         for(int i = 0, idx = 0; i < _octave_num; i++)
01066         {
01067                 CuTexImage * got = GetBaseLevel(i + _octave_min, DATA_GRAD) +  1;
01068                 //compute the gradient
01069                 for(int j = 0; j <  param._dog_level_num ; j++, got++, idx++)
01070                 {
01071                         if(_levelFeatureNum[idx] > 0)   got->CopyToTexture2D();
01072                 }
01073         }
01074         if(GlobalUtil::_timingS)
01075         {
01076                 ProgramCU::FinishCUDA();
01077                 if(GlobalUtil::_verbose)
01078                 {
01079                         t1 = CLOCK();
01080                         std::cout       <<"<Copy Grad/Orientation>\t"<<(t1-ts)<<"\n";
01081                 }
01082         }
01083 }
01084 
01085 void PyramidCU::ComputeGradient() 
01086 {
01087 
01088         int i, j;
01089         double ts, t1;
01090 
01091         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
01092 
01093         for(i = _octave_min; i < _octave_min + _octave_num; i++)
01094         {
01095                 CuTexImage * gus = GetBaseLevel(i) +  1;
01096                 CuTexImage * dog = GetBaseLevel(i, DATA_DOG) +  1;
01097                 CuTexImage * got = GetBaseLevel(i, DATA_GRAD) +  1;
01098 
01099                 //compute the gradient
01100                 for(j = 0; j <  param._dog_level_num ; j++, gus++, dog++, got++)
01101                 {
01102                         ProgramCU::ComputeDOG(gus, dog, got);
01103                 }
01104         }
01105         if(GlobalUtil::_timingS)
01106         {
01107                 ProgramCU::FinishCUDA();
01108                 if(GlobalUtil::_verbose)
01109                 {
01110                         t1 = CLOCK();
01111                         std::cout       <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n";
01112                 }
01113         }
01114 }
01115 
01116 int PyramidCU::FitHistogramPyramid(CuTexImage* tex)
01117 {
01118         CuTexImage *htex;
01119         int hist_level_num = _hpLevelNum - _pyramid_octave_first / 2; 
01120         htex = _histoPyramidTex + hist_level_num - 1;
01121         int w = (tex->GetImgWidth() + 2) >> 2;
01122         int h = tex->GetImgHeight();
01123         int count = 0; 
01124         for(int k = 0; k < hist_level_num; k++, htex--)
01125         {
01126                 //htex->SetImageSize(w, h);     
01127                 htex->InitTexture(w, h, 4); 
01128                 ++count;
01129                 if(w == 1)
01130             break;
01131                 w = (w + 3)>>2; 
01132         }
01133         return count;
01134 }
01135 
01136 void PyramidCU::GetFeatureOrientations() 
01137 {
01138 
01139         CuTexImage * ftex = _featureTex;
01140         int * count      = _levelFeatureNum;
01141         float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num);
01142 
01143         for(int i = 0; i < _octave_num; i++)
01144         {
01145                 CuTexImage* got = GetBaseLevel(i + _octave_min, DATA_GRAD) + 1;
01146                 CuTexImage* key = GetBaseLevel(i + _octave_min, DATA_KEYPOINT) + 2;
01147 
01148                 for(int j = 0; j < param._dog_level_num; j++, ftex++, count++, got++, key++)
01149                 {
01150                         if(*count<=0)continue;
01151 
01152                         //if(ftex->GetImgWidth() < *count) ftex->InitTexture(*count, 1, 4);
01153 
01154                         sigma = param.GetLevelSigma(j+param._level_min+1);
01155 
01156                         ProgramCU::ComputeOrientation(ftex, got, key, sigma, sigma_step, _existing_keypoints);          
01157                 }
01158         }
01159 
01160         if(GlobalUtil::_timingS)ProgramCU::FinishCUDA();
01161         if(ProgramCU::CheckErrorCUDA("PyramidCU::GetFeatureOrientations")) SetFailStatus();
01162 
01163 }
01164 
01165 void PyramidCU::GetSimplifiedOrientation() 
01166 {
01167         //no simplified orientation
01168         GetFeatureOrientations();
01169 }
01170 
01171 CuTexImage* PyramidCU::GetBaseLevel(int octave, int dataName)
01172 {
01173         if(octave <_octave_min || octave > _octave_min + _octave_num) return NULL;
01174         int offset = (_pyramid_octave_first + octave - _octave_min) * param._level_num;
01175         int num = param._level_num * _pyramid_octave_num;
01176         if (dataName == DATA_ROT) dataName = DATA_GRAD;
01177         return _allPyramid + num * dataName + offset;
01178 }
01179 
01180 #endif
01181 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines


rgbd_registration
Author(s): Ross Kidson
autogenerated on Thu May 23 2013 15:36:53