PyramidGL.cpp
Go to the documentation of this file.
00001 
00002 //      File:           PyramidGL.cpp
00003 //      Author:         Changchang Wu
00004 //      Description : implementation of PyramidGL/PyramidNaive/PyramidPackdc .
00005 //
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 #include "GL/glew.h"
00024 #include <iostream>
00025 #include <iomanip>
00026 #include <vector>
00027 #include <algorithm>
00028 #include <fstream>
00029 #include <math.h>
00030 #include <string.h>
00031 using namespace std;
00032 
00033 #include "GlobalUtil.h"
00034 #include "GLTexImage.h"
00035 #include "SiftGPU.h"
00036 #include "ShaderMan.h"
00037 #include "SiftPyramid.h"
00038 #include "ProgramGLSL.h"
00039 #include "PyramidGL.h"
00040 #include "FrameBufferObject.h"
00041 
00042 
00043 #if defined(__SSE__) || _MSC_VER > 1200
00044 #define USE_SSE_FOR_SIFTGPU
00045 #include <xmmintrin.h>
00046 #else
00047 //#pragma message( "SSE optimization off!\n" )
00048 #endif
00049 
00050 
00051 #define USE_TIMING()            double t, t0, tt;
00052 #define OCTAVE_START()          if(GlobalUtil::_timingO){       t = t0 = CLOCK();       cout<<"#"<<i+_down_sample_factor<<"\t"; }
00053 #define LEVEL_FINISH()          if(GlobalUtil::_timingL){       glFinish();     tt = CLOCK();cout<<(tt-t)<<"\t";        t = CLOCK();}
00054 #define OCTAVE_FINISH()         if(GlobalUtil::_timingO)cout<<"|\t"<<(CLOCK()-t0)<<endl;
00055 
00056 
00058 // Construction/Destruction
00060 PyramidNaive::PyramidNaive(SiftParam& sp): PyramidGL(sp)
00061 {
00062         _texPyramid = NULL;
00063         _auxPyramid = NULL;
00064 }
00065 
00066 PyramidNaive::~PyramidNaive()
00067 {
00068         DestroyPyramidData();
00069 }
00070 
00071 //align must be 2^i
00072 void PyramidGL::	GetAlignedStorageSize(int num, int align,  int &fw, int &fh)
00073 {
00074         if(num <=0)
00075         {
00076                 fw = fh = 0;
00077         }else if(num < align*align)
00078         {
00079                 fw = align;
00080                 fh = (int)ceil(double(num) / fw);
00081         }else if(GlobalUtil::_NarrowFeatureTex)
00082         {
00083                 double dn = double(num);
00084                 int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/align);   
00085                 fw = align * nb;        
00086                 fh = (int)ceil(dn /fw);
00087         }else
00088         {
00089                 double dn = double(num);
00090                 int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/align);
00091                 fh = align * nb;
00092                 if(nb <=1)
00093                 {
00094                         fw = (int)ceil(dn / fh);
00095                         //align this dimension to blocksize
00096                         fw = ((int) ceil(double(fw) /align))*align;
00097                 }else
00098                 {
00099                         fw = GlobalUtil::_texMaxDim;
00100                 }
00101 
00102         }
00103 
00104 
00105 }
00106 
00107 void PyramidGL::GetTextureStorageSize(int num, int &fw, int& fh)
00108 {
00109         if(num <=0)
00110         {
00111                 fw = fh = 0;
00112         }else if(num <= GlobalUtil::_FeatureTexBlock)
00113         {
00114                 fw = num;
00115                 fh = 1;
00116         }else if(GlobalUtil::_NarrowFeatureTex)
00117         {
00118                 double dn = double(num);
00119                 int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/GlobalUtil::_FeatureTexBlock);    
00120                 fw = GlobalUtil::_FeatureTexBlock * nb; 
00121                 fh = (int)ceil(dn /fw);
00122         }else
00123         {
00124                 double dn = double(num);
00125                 int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/GlobalUtil::_FeatureTexBlock);
00126                 fh = GlobalUtil::_FeatureTexBlock * nb;
00127                 if(nb <=1)
00128                 {
00129                         fw = (int)ceil(dn / fh);
00130 
00131                         //align this dimension to blocksize
00132 
00133                         //
00134                         if( fw < fh)
00135                         {
00136                                 int temp = fh;
00137                                 fh = fw;
00138                                 fw = temp;
00139                         }
00140                 }else
00141                 {
00142                         fw = GlobalUtil::_texMaxDim;
00143                 }
00144         }
00145 }
00146 
00147 void PyramidNaive::DestroyPyramidData()
00148 {
00149         if(_texPyramid)
00150         {
00151                 delete [] _texPyramid;
00152                 _texPyramid = NULL;
00153         }
00154         if(_auxPyramid)
00155         {
00156                 delete [] _auxPyramid;  
00157                 _auxPyramid = NULL;
00158         }
00159 }
00160 PyramidGL::PyramidGL(SiftParam &sp):SiftPyramid(sp)
00161 {
00162         _featureTex = NULL;
00163         _orientationTex = NULL;
00164         _descriptorTex = NULL;
00165         _histoPyramidTex = NULL;        
00167 
00168     InitializeContext();
00169 }
00170 
00171 PyramidGL::~PyramidGL()
00172 {
00173         DestroyPerLevelData();
00174         DestroySharedData();
00175         ShaderMan::DestroyShaders();
00176 }
00177 
00178 void PyramidGL::InitializeContext()
00179 {
00180     GlobalUtil::InitGLParam(0);
00181     if(!GlobalUtil::_GoodOpenGL) return;
00182 
00184         ShaderMan::InitShaderMan(param);
00185 }
00186 
00187 void PyramidGL::DestroyPerLevelData()
00188 {
00189         //integers vector to store the feature numbers.
00190         if(_levelFeatureNum)
00191         {
00192                 delete [] _levelFeatureNum;
00193                 _levelFeatureNum = NULL;
00194         }
00195         //texture used to store features
00196         if(     _featureTex)
00197         {
00198                 delete [] _featureTex;
00199                 _featureTex =   NULL;
00200         }
00201         //texture used for multi-orientation 
00202         if(_orientationTex)
00203         {
00204                 delete [] _orientationTex;
00205                 _orientationTex = NULL;
00206         }
00207         int no = _octave_num* param._dog_level_num;
00208 
00209         //two sets of vbos used to display the features
00210         if(_featureDisplayVBO)
00211         {
00212                 glDeleteBuffers(no, _featureDisplayVBO);
00213                 delete [] _featureDisplayVBO;
00214                 _featureDisplayVBO = NULL;
00215         }
00216         if( _featurePointVBO)
00217         {
00218                 glDeleteBuffers(no, _featurePointVBO);
00219                 delete [] _featurePointVBO;
00220                 _featurePointVBO = NULL;
00221         }
00222 
00223 }
00224 
00225 void PyramidGL::DestroySharedData()
00226 {
00227         //histogram reduction
00228         if(_histoPyramidTex)
00229         {
00230                 delete[]        _histoPyramidTex;
00231                 _hpLevelNum = 0;
00232                 _histoPyramidTex = NULL;
00233         }
00234 
00235         //descriptor storage shared by all levels
00236         if(_descriptorTex)
00237         {
00238                 delete [] _descriptorTex;
00239                 _descriptorTex = NULL;
00240         }
00241         //cpu reduction buffer.
00242         if(_histo_buffer)
00243         {
00244                 delete[] _histo_buffer;
00245                 _histo_buffer = 0;
00246         }
00247 }
00248 
00249 void PyramidNaive::FitHistogramPyramid()
00250 {
00251         GLTexImage * tex, *htex;
00252         int hist_level_num = _hpLevelNum - _pyramid_octave_first; 
00253 
00254         tex = GetBaseLevel(_octave_min , DATA_KEYPOINT) + 2;
00255         htex = _histoPyramidTex + hist_level_num - 1;
00256         int w = tex->GetImgWidth() >> 1;
00257         int h = tex->GetImgHeight() >> 1;
00258 
00259         for(int k = 0; k <hist_level_num -1; k++, htex--)
00260         {
00261                 if(htex->GetImgHeight()!= h || htex->GetImgWidth() != w)
00262                 {       
00263                         htex->SetImageSize(w, h);
00264                         htex->ZeroHistoMargin();
00265                 }
00266 
00267                 w = (w + 1)>>1; h = (h + 1) >> 1;
00268         }
00269 }
00270 
00271 void PyramidNaive::FitPyramid(int w, int h)
00272 {
00273         //(w, h) <= (_pyramid_width, _pyramid_height);
00274 
00275         _pyramid_octave_first = 0;
00276         //
00277         _octave_num  = GlobalUtil::_octave_num_default;
00278 
00279         int _octave_num_max = GetRequiredOctaveNum(min(w, h));
00280 
00281         if(_octave_num < 1 || _octave_num > _octave_num_max) 
00282         {
00283                 _octave_num = _octave_num_max;
00284         }
00285 
00286 
00287         int pw = _pyramid_width>>1, ph = _pyramid_height>>1;
00288         while(_pyramid_octave_first + _octave_num < _pyramid_octave_num &&  
00289                 pw >= w && ph >= h)
00290         {
00291                 _pyramid_octave_first++;
00292                 pw >>= 1;
00293                 ph >>= 1;
00294         }
00295 
00296         for(int i = 0; i < _octave_num; i++)
00297         {
00298                 GLTexImage * tex = GetBaseLevel(i + _octave_min);
00299                 GLTexImage * aux = GetBaseLevel(i + _octave_min, DATA_KEYPOINT);
00300                 for(int j = param._level_min; j <= param._level_max; j++, tex++, aux++)
00301                 {
00302                         tex->SetImageSize(w, h);
00303                         aux->SetImageSize(w, h);
00304                 }
00305                 w>>=1;
00306                 h>>=1;
00307         }
00308 }
00309 void PyramidNaive::InitPyramid(int w, int h, int ds)
00310 {
00311         int wp, hp, toobig = 0;
00312         if(ds == 0)
00313         {
00314                 _down_sample_factor = 0;
00315                 if(GlobalUtil::_octave_min_default>=0)
00316                 {
00317                         wp = w >> GlobalUtil::_octave_min_default;
00318                         hp = h >> GlobalUtil::_octave_min_default;
00319                 }else 
00320                 {
00321                         wp = w << (-GlobalUtil::_octave_min_default);
00322                         hp = h << (-GlobalUtil::_octave_min_default);
00323                 }
00324                 _octave_min = _octave_min_default;
00325         }else
00326         {
00327                 //must use 0 as _octave_min; 
00328                 _octave_min = 0;
00329                 _down_sample_factor = ds;
00330                 w >>= ds;
00331                 h >>= ds;
00332                 wp = w;
00333                 hp = h; 
00334 
00335         }
00336 
00337         while(wp > GlobalUtil::_texMaxDim || hp > GlobalUtil::_texMaxDim)
00338         {
00339                 _octave_min ++;
00340                 wp >>= 1;
00341                 hp >>= 1;
00342                 toobig = 1;
00343         }
00344 
00345         if(toobig && GlobalUtil::_verbose && _octave_min > 0)
00346         {
00347 
00348                 std::cout<< "**************************************************************\n"
00349                                         "Image larger than allowed dimension, data will be downsampled!\n"
00350                                         "use -maxd to change the settings\n"
00351                                         "***************************************************************\n";
00352         }
00353 
00354         if( wp == _pyramid_width && hp == _pyramid_height && _allocated )
00355         {
00356                 FitPyramid(wp, hp);
00357         }else if(GlobalUtil::_ForceTightPyramid || _allocated ==0)
00358         {
00359                 ResizePyramid(wp, hp);
00360         }
00361         else if( wp > _pyramid_width || hp > _pyramid_height )
00362         {
00363                 ResizePyramid(max(wp, _pyramid_width), max(hp, _pyramid_height));
00364                 if(wp < _pyramid_width || hp < _pyramid_height)  FitPyramid(wp, hp);
00365         }
00366         else
00367         {
00368                 //try use the pyramid allocated for large image on small input images
00369                 FitPyramid(wp, hp);
00370         }
00371 
00372         //select the initial smoothing filter according to the new _octave_min
00373         ShaderMan::SelectInitialSmoothingFilter(_octave_min + _down_sample_factor, param);
00374 }
00375 
00376 void PyramidNaive::ResizePyramid( int w,  int h)
00377 {
00378         //
00379         unsigned int totalkb = 0;
00380         int _octave_num_new, input_sz;
00381         int i, j;
00382         GLTexImage * tex, *aux;
00383         //
00384 
00385         if(_pyramid_width == w && _pyramid_height == h && _allocated) return;
00386 
00387         if(w > GlobalUtil::_texMaxDim || h > GlobalUtil::_texMaxDim) return ;
00388 
00389         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<w<<"x"<<h<<endl;
00390         //first octave does not change
00391         _pyramid_octave_first = 0;
00392 
00393         
00394         //compute # of octaves
00395 
00396         input_sz = min(w,h) ;
00397 
00398 
00399         _pyramid_width =  w;
00400         _pyramid_height =  h;
00401 
00402         //reset to preset parameters
00403         _octave_num_new  = GlobalUtil::_octave_num_default;
00404 
00405         if(_octave_num_new < 1) _octave_num_new = GetRequiredOctaveNum(input_sz)  ;
00406 
00407         if(_pyramid_octave_num != _octave_num_new)
00408         {
00409                 //destroy the original pyramid if the # of octave changes
00410                 if(_octave_num >0)
00411                 {
00412                         DestroyPerLevelData();
00413                         DestroyPyramidData();
00414                 }
00415                 _pyramid_octave_num = _octave_num_new;
00416         }
00417 
00418         _octave_num = _pyramid_octave_num;
00419 
00420         int noct = _octave_num;
00421         int nlev = param._level_num;
00422 
00423         //      //initialize the pyramid
00424         if(_texPyramid==NULL)   _texPyramid = new GLTexImage[ noct* nlev ];
00425         if(_auxPyramid==NULL)   _auxPyramid = new GLTexImage[ noct* nlev ];
00426 
00427 
00428         tex = GetBaseLevel(_octave_min, DATA_GAUSSIAN);
00429         aux = GetBaseLevel(_octave_min, DATA_KEYPOINT);
00430         for(i = 0; i< noct; i++)
00431         {
00432                 totalkb += (nlev * w * h * 16 / 1024);
00433                 for( j = 0; j< nlev; j++, tex++)
00434                 {
00435                         tex->InitTexture(w, h);
00436                         //tex->AttachToFBO(0);
00437                 }
00438                 //several auxilary textures are not actually required
00439                 totalkb += ((nlev - 3) * w * h * 16 /1024);
00440                 for( j = 0; j< nlev ; j++, aux++)
00441                 {
00442                         if(j < 2) continue;
00443                         if(j >= nlev - 1) continue;
00444                         aux->InitTexture(w, h, 0);
00445                         //aux->AttachToFBO(0);
00446                 }
00447 
00448                 w>>=1;
00449                 h>>=1;
00450         }
00451 
00452         totalkb += ResizeFeatureStorage();
00453 
00454 
00455         //
00456         _allocated = 1;
00457 
00458         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<(totalkb/1024)<<"MB\n";
00459 
00460 }
00461 
00462 
00463 int PyramidGL::ResizeFeatureStorage()
00464 {
00465         int totalkb = 0;
00466         if(_levelFeatureNum==NULL)      _levelFeatureNum = new int[_octave_num * param._dog_level_num];
00467         std::fill(_levelFeatureNum, _levelFeatureNum+_octave_num * param._dog_level_num, 0); 
00468 
00469         int wmax = GetBaseLevel(_octave_min)->GetDrawWidth();
00470         int hmax = GetBaseLevel(_octave_min)->GetDrawHeight();
00471         int w ,h, i;
00472 
00473         //use a fbo to initialize textures..
00474         FrameBufferObject fbo;
00475         
00476         //
00477         if(_histo_buffer == NULL) _histo_buffer = new float[1 << (2 + 2 * GlobalUtil::_ListGenSkipGPU)];
00478         //histogram for feature detection
00479 
00480         int num = (int)ceil(log(double(max(wmax, hmax)))/log(2.0));
00481 
00482         if( _hpLevelNum != num)
00483         {
00484                 _hpLevelNum = num;
00485                 if(GlobalUtil::_ListGenGPU)
00486                 {
00487                         if(_histoPyramidTex ) delete [] _histoPyramidTex;
00488                         _histoPyramidTex = new GLTexImage[_hpLevelNum];
00489                         w = h = 1 ;
00490                         for(i = 0; i < _hpLevelNum; i++)
00491                         {
00492                                 _histoPyramidTex[i].InitTexture(w, h, 0);
00493                                 _histoPyramidTex[i].AttachToFBO(0);
00494                                 w<<=1;
00495                                 h<<=1;
00496                         }
00497                 }
00498         }
00499 
00500         // (4 ^ (_hpLevelNum) -1 / 3) pixels
00501         if(GlobalUtil::_ListGenGPU) totalkb += (((1 << (2 * _hpLevelNum)) -1) / 3 * 16 / 1024);
00502 
00503 
00504 
00505         //initialize the feature texture
00506 
00507         int idx = 0, n = _octave_num * param._dog_level_num;
00508         if(_featureTex==NULL)   _featureTex = new GLTexImage[n];
00509         if(GlobalUtil::_MaxOrientation >1 && GlobalUtil::_OrientationPack2==0)
00510         {
00511                 if(_orientationTex== NULL)              _orientationTex = new GLTexImage[n];
00512         }
00513 
00514 
00515         for(i = 0; i < _octave_num; i++)
00516         {
00517                 GLTexImage * tex = GetBaseLevel(i+_octave_min);
00518                 int fmax = int(tex->GetImgWidth()*tex->GetImgHeight()*GlobalUtil::_MaxFeaturePercent);
00519                 int fw, fh;
00520                 //
00521                 if(fmax > GlobalUtil::_MaxLevelFeatureNum) fmax = GlobalUtil::_MaxLevelFeatureNum;
00522                 else if(fmax < 32) fmax = 32;   //give it at least a space of 32 feature
00523 
00524                 GetTextureStorageSize(fmax, fw, fh);
00525                 
00526                 for(int j = 0; j < param._dog_level_num; j++, idx++)
00527                 {
00528 
00529                         _featureTex[idx].InitTexture(fw, fh, 0);
00530                         _featureTex[idx].AttachToFBO(0);
00531                         //
00532                         if(_orientationTex)
00533                         {
00534                                 _orientationTex[idx].InitTexture(fw, fh, 0);
00535                                 _orientationTex[idx].AttachToFBO(0);
00536                         }
00537                 }
00538                 totalkb += fw * fh * 16 * param._dog_level_num * (_orientationTex? 2 : 1) /1024;
00539         }
00540 
00541 
00542         //this just need be initialized once
00543         if(_descriptorTex==NULL)
00544         {
00545                 //initialize feature texture pyramid
00546                 wmax = _featureTex->GetImgWidth();
00547                 hmax = _featureTex->GetImgHeight();
00548 
00549                 int nf, ns;
00550                 if(GlobalUtil::_DescriptorPPT)
00551                 {
00552                         //32*4 = 128. 
00553                         nf = 32 / GlobalUtil::_DescriptorPPT;   // how many textures we need
00554                         ns = max(4, GlobalUtil::_DescriptorPPT);                    // how many point in one texture for one descriptor
00555                 }else
00556                 {
00557                         //at least one, resue for visualization and other work
00558                         nf = 1; ns = 4;
00559                 }
00560                 //
00561                 _alignment = ns;
00562                 //
00563                 _descriptorTex = new GLTexImage[nf];
00564 
00565                 int fw, fh;
00566                 GetAlignedStorageSize(hmax*wmax* max(ns, 10), _alignment, fw, fh);
00567 
00568                 if(fh < hmax ) fh = hmax;
00569                 if(fw < wmax ) fw = wmax;
00570 
00571                 totalkb += ( fw * fh * nf * 16 /1024);
00572                 for(i =0; i < nf; i++)
00573                 {
00574                         _descriptorTex[i].InitTexture(fw, fh);
00575                 }
00576         }else
00577         {
00578                 int nf = GlobalUtil::_DescriptorPPT? 32 / GlobalUtil::_DescriptorPPT: 1;
00579                 totalkb += nf * _descriptorTex[0].GetTexWidth() * _descriptorTex[0].GetTexHeight() * 16 /1024;
00580         }
00581         return totalkb;
00582 }
00583 
00584 
00585 void PyramidNaive::BuildPyramid(GLTexInput *input)
00586 {
00587         USE_TIMING();
00588         GLTexPacked * tex;
00589         FilterProgram** filter;
00590         FrameBufferObject fbo;
00591 
00592         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
00593         input->FitTexViewPort();
00594 
00595         for (int i = _octave_min; i < _octave_min + _octave_num; i++)
00596         {
00597 
00598                 tex = (GLTexPacked*)GetBaseLevel(i);
00599                 filter = ShaderMan::s_bag->f_gaussian_step;
00600 
00601                 OCTAVE_START();
00602 
00603                 if( i == _octave_min )
00604                 {
00605                         if(i < 0)       TextureUpSample(tex, input, 1<<(-i)     );                      
00606                         else        TextureDownSample(tex, input, 1<<i);
00607                 ShaderMan::FilterInitialImage(tex, NULL);
00608                 }else
00609                 {
00610                         TextureDownSample(tex, GetLevelTexture(i-1, param._level_ds)); 
00611             ShaderMan::FilterSampledImage(tex, NULL); 
00612         }
00613                 LEVEL_FINISH();
00614 
00615                 for(int j = param._level_min + 1; j <=  param._level_max ; j++, tex++, filter++)
00616                 {
00617                         // filtering
00618             ShaderMan::FilterImage(*filter, tex+1, tex, NULL);
00619                         LEVEL_FINISH();
00620                 }
00621                 OCTAVE_FINISH();
00622 
00623         }
00624         if(GlobalUtil::_timingS)        glFinish();
00625         UnloadProgram();
00626 }
00627 
00628 
00629 
00630 
00631 
00632 
00633 GLTexImage*  PyramidNaive::GetLevelTexture(int octave, int level, int dataName)
00634 {
00635         if(octave <_octave_min || octave > _octave_min + _octave_num) return NULL;
00636         switch(dataName)
00637         {
00638                 case DATA_GAUSSIAN:
00639                 case DATA_DOG:
00640                 case DATA_GRAD:
00641                 case DATA_ROT:
00642                         return _texPyramid+ (_pyramid_octave_first + octave - _octave_min) * param._level_num + (level - param._level_min);
00643                 case DATA_KEYPOINT:
00644                         return _auxPyramid + (_pyramid_octave_first + octave - _octave_min) * param._level_num + (level - param._level_min);
00645                 default:
00646                         return NULL;
00647         }
00648 }
00649 
00650 GLTexImage*  PyramidNaive::GetLevelTexture(int octave, int level)
00651 {
00652         return _texPyramid+ (_pyramid_octave_first + octave - _octave_min) * param._level_num 
00653                 + (level - param._level_min);
00654 }
00655 
00656 //in the packed implementation
00657 // DATA_GAUSSIAN, DATA_DOG, DATA_GAD will be stored in different textures.
00658 GLTexImage*  PyramidNaive::GetBaseLevel(int octave, int dataName)
00659 {
00660         if(octave <_octave_min || octave > _octave_min + _octave_num) return NULL;
00661         switch(dataName)
00662         {
00663                 case DATA_GAUSSIAN:
00664                 case DATA_DOG:
00665                 case DATA_GRAD:
00666                 case DATA_ROT:
00667                         return _texPyramid+ (_pyramid_octave_first + octave - _octave_min) * param._level_num;
00668                 case DATA_KEYPOINT:
00669                         return _auxPyramid + (_pyramid_octave_first + octave - _octave_min) * param._level_num;
00670                 default:
00671                         return NULL;
00672         }
00673 }
00674 
00675 
00676 
00677 
00678 
00679 
00680 
00681 
00682 
00683 void PyramidNaive::ComputeGradient()
00684 {
00685 
00686         int i, j;
00687         double  ts, t1;
00688         GLTexImage * tex;
00689         FrameBufferObject fbo;
00690 
00691 
00692         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
00693         
00694         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
00695 
00696         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
00697         {
00698                 for( j = param._level_min + 1 ; j < param._level_max ; j++)
00699                 {
00700                         tex = GetLevelTexture(i, j);
00701                         tex->FitTexViewPort();
00702                         tex->AttachToFBO(0);
00703                         tex->BindTex();
00704                         ShaderMan::UseShaderGradientPass();
00705                         tex->DrawQuadMT4();
00706                 }
00707         }               
00708 
00709         if(GlobalUtil::_timingS && GlobalUtil::_verbose)
00710         {
00711                 glFinish();
00712                 t1 = CLOCK();   
00713                 std::cout<<"<Compute Gradient>\t"<<(t1-ts)<<"\n";
00714         }
00715 
00716         UnloadProgram();
00717         GLTexImage::UnbindMultiTex(3);
00718         fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
00719 }
00720 
00721 
00722 //keypoint detection with subpixel localization
00723 void PyramidNaive::DetectKeypointsEX()
00724 {
00725         int i, j;
00726         double t0, t, ts, t1, t2;
00727         GLTexImage * tex, *aux;
00728         FrameBufferObject fbo;
00729 
00730         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
00731         
00732         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
00733         //extra gradient data required for visualization
00734         int gradient_only_levels[2] = {param._level_min +1, param._level_max};
00735         int n_gradient_only_level = GlobalUtil::_UseSiftGPUEX ? 2 : 1;
00736         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
00737         {
00738                 for( j =0; j < n_gradient_only_level ; j++)
00739                 {
00740                         tex = GetLevelTexture(i, gradient_only_levels[j]);
00741                         tex->FitTexViewPort();
00742                         tex->AttachToFBO(0);
00743                         tex->BindTex();
00744                         ShaderMan::UseShaderGradientPass();
00745                         tex->DrawQuadMT4();
00746                 }
00747         }               
00748 
00749         if(GlobalUtil::_timingS && GlobalUtil::_verbose)
00750         {
00751                 glFinish();
00752                 t1 = CLOCK();
00753         }
00754 
00755         GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT };
00756         glDrawBuffers(2, buffers);
00757         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
00758         {
00759                 if(GlobalUtil::_timingO)
00760                 {
00761                         t0 = CLOCK();
00762                         std::cout<<"#"<<(i + _down_sample_factor)<<"\t";
00763                 }
00764                 tex = GetBaseLevel(i) + 2;
00765                 aux = GetBaseLevel(i, DATA_KEYPOINT) +2;
00766                 aux->FitTexViewPort();
00767 
00768                 for( j = param._level_min + 2; j <  param._level_max ; j++, aux++, tex++)
00769                 {
00770                         if(GlobalUtil::_timingL)t = CLOCK();            
00771                         tex->AttachToFBO(0);
00772                         aux->AttachToFBO(1);
00773                         glActiveTexture(GL_TEXTURE0);
00774                         tex->BindTex();
00775                         glActiveTexture(GL_TEXTURE1);
00776                         (tex+1)->BindTex();
00777                         glActiveTexture(GL_TEXTURE2);
00778                         (tex-1)->BindTex();
00779                         ShaderMan::UseShaderKeypoint((tex+1)->GetTexID(), (tex-1)->GetTexID());
00780                         aux->DrawQuadMT8();
00781         
00782                         if(GlobalUtil::_timingL)
00783                         {
00784                                 glFinish();
00785                                 std::cout<<(CLOCK()-t)<<"\t";
00786                         }
00787                         tex->DetachFBO(0);
00788                         aux->DetachFBO(1);
00789                 }
00790                 if(GlobalUtil::_timingO)
00791                 {
00792                         std::cout<<"|\t"<<(CLOCK()-t0)<<"\n";
00793                 }
00794         }
00795 
00796         if(GlobalUtil::_timingS)
00797         {
00798                 glFinish();
00799                 t2 = CLOCK();
00800                 if(GlobalUtil::_verbose) 
00801                         std::cout       <<"<Get Keypoints ..  >\t"<<(t2-t1)<<"\n"
00802                                                 <<"<Extra Gradient..  >\t"<<(t1-ts)<<"\n";
00803         }
00804         UnloadProgram();
00805         GLTexImage::UnbindMultiTex(3);
00806         fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
00807 
00808 
00809 }
00810 
00811 void PyramidNaive::GenerateFeatureList(int i, int j)
00812 {
00813         int hist_level_num = _hpLevelNum - _pyramid_octave_first; 
00814         int hist_skip_gpu = GlobalUtil::_ListGenSkipGPU; 
00815     int idx = i * param._dog_level_num + j;
00816     GLTexImage* htex, *ftex, *tex;
00817         tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2 + j;
00818     ftex = _featureTex + idx;
00819         htex = _histoPyramidTex + hist_level_num - 1 - i;
00820 
00822         glActiveTexture(GL_TEXTURE0);
00823         tex->BindTex();
00824         htex->AttachToFBO(0);
00825         int tight = ((htex->GetImgWidth() * 2 == tex->GetImgWidth() -1 || tex->GetTexWidth() == tex->GetImgWidth()) &&
00826                                  (htex->GetImgHeight() *2 == tex->GetImgHeight()-1 || tex->GetTexHeight() == tex->GetImgHeight()));
00827         ShaderMan::UseShaderGenListInit(tex->GetImgWidth(), tex->GetImgHeight(), tight);
00828         htex->FitTexViewPort();
00829         //this uses the fact that no feature is on the edge.
00830         htex->DrawQuadReduction();
00831 
00832         //reduction..
00833         htex--;
00834 
00835         //this part might have problems on several GPUS
00836         //because the output of one pass is the input of the next pass
00837         //need to call glFinish to make it right
00838         //but too much glFinish makes it slow
00839         for(int k = 0; k <hist_level_num - i - 1 - hist_skip_gpu; k++, htex--)
00840         {
00841                 htex->AttachToFBO(0);
00842                 htex->FitTexViewPort();
00843                 (htex+1)->BindTex();
00844                 ShaderMan::UseShaderGenListHisto();
00845                 htex->DrawQuadReduction();                                      
00846         }
00847 
00848         //
00849         if(hist_skip_gpu == 0)
00850         {       
00851                 //read back one pixel
00852                 float fn[4], fcount;
00853                 glReadPixels(0, 0, 1, 1, GL_RGBA , GL_FLOAT, fn);
00854                 fcount = (fn[0] + fn[1] + fn[2] + fn[3]);
00855                 if(fcount < 1) fcount = 0;
00856 
00857 
00858                 _levelFeatureNum[ idx] = (int)(fcount);
00859                 SetLevelFeatureNum(idx, (int)fcount);
00860                 _featureNum += int(fcount);
00861 
00862                 //
00863                 if(fcount < 1.0) return;
00864         
00865 
00867 
00868                 htex=  _histoPyramidTex;
00869 
00870                 htex->BindTex();
00871 
00872                 //first pass
00873                 ftex->AttachToFBO(0);
00874                 if(GlobalUtil::_MaxOrientation>1)
00875                 {
00876                         //this is very important...
00877                         ftex->FitRealTexViewPort();
00878                         glClear(GL_COLOR_BUFFER_BIT);
00879                         glFinish();
00880                 }else
00881                 {
00882                         ftex->FitTexViewPort();
00883             //glFinish();
00884                 }
00885 
00886 
00887                 ShaderMan::UseShaderGenListStart((float)ftex->GetImgWidth(), htex->GetTexID());
00888 
00889                 ftex->DrawQuad();
00890                 //make sure it finishes before the next step
00891                 ftex->DetachFBO(0);
00892 
00893                 //pass on each pyramid level
00894                 htex++;
00895         }else
00896         {
00897 
00898                 int tw = htex[1].GetDrawWidth(), th = htex[1].GetDrawHeight();
00899                 int fc = 0;
00900                 glReadPixels(0, 0, tw, th, GL_RGBA , GL_FLOAT, _histo_buffer);  
00901                 _keypoint_buffer.resize(0);
00902                 for(int y = 0, pos = 0; y < th; y++)
00903                 {
00904                         for(int x= 0; x < tw; x++)
00905                         {
00906                                 for(int c = 0; c < 4; c++, pos++)
00907                                 {
00908                                         int ss =  (int) _histo_buffer[pos]; 
00909                                         if(ss == 0) continue;
00910                                         float ft[4] = {2 * x + (c%2? 1.5f:  0.5f), 2 * y + (c>=2? 1.5f: 0.5f), 0, 1 };
00911                                         for(int t = 0; t < ss; t++)
00912                                         {
00913                                                 ft[2] = (float) t; 
00914                                                 _keypoint_buffer.insert(_keypoint_buffer.end(), ft, ft+4);
00915                                         }
00916                                         fc += (int)ss; 
00917                                 }
00918                         }
00919                 }
00920                 _levelFeatureNum[ idx] = fc;
00921                 SetLevelFeatureNum(idx, fc);
00922                 if(fc == 0)  return;
00923         _featureNum += fc;
00925                 ftex->AttachToFBO(0);
00926                 if(GlobalUtil::_MaxOrientation>1)
00927                 {
00928                         ftex->FitRealTexViewPort();
00929                         glClear(GL_COLOR_BUFFER_BIT);
00930                 }else
00931                 {                                       
00932                         ftex->FitTexViewPort();
00933                 }
00934                 _keypoint_buffer.resize(ftex->GetDrawWidth() * ftex->GetDrawHeight()*4, 0);
00936                 glActiveTexture(GL_TEXTURE0);
00937                 ftex->BindTex();
00938                 glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, ftex->GetDrawWidth(),
00939                         ftex->GetDrawHeight(), GL_RGBA, GL_FLOAT, &_keypoint_buffer[0]);
00940                 htex += 2;
00941         }
00942 
00943         for(int lev = 1 + hist_skip_gpu; lev < hist_level_num  - i; lev++, htex++)
00944         {
00945 
00946                 glActiveTexture(GL_TEXTURE0);
00947                 ftex->BindTex();
00948                 ftex->AttachToFBO(0);
00949                 glActiveTexture(GL_TEXTURE1);
00950                 htex->BindTex();
00951                 ShaderMan::UseShaderGenListStep(ftex->GetTexID(), htex->GetTexID());
00952                 ftex->DrawQuad();
00953                 ftex->DetachFBO(0);     
00954         }
00955         GLTexImage::UnbindMultiTex(2);
00956 
00957 }
00958 
00959 //generate feature list on GPU
00960 void PyramidNaive::GenerateFeatureList()
00961 {
00962         //generate the histogram0pyramid
00963         FrameBufferObject fbo;
00964         glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
00965         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
00966         double t1, t2; 
00967         int ocount, reverse = (GlobalUtil::_TruncateMethod == 1);
00968         _featureNum = 0;
00969 
00970         FitHistogramPyramid();
00971 
00972         //for(int i = 0, idx = 0; i < _octave_num; i++)
00973     FOR_EACH_OCTAVE(i, reverse)
00974         {
00975                 //output
00976                 if(GlobalUtil::_timingO)
00977                 {
00978             t1= CLOCK();
00979                         ocount = 0;
00980                         std::cout<<"#"<<i+_octave_min + _down_sample_factor<<":\t";
00981                 }
00982                 //for(int j = 0; j < param._dog_level_num; j++, idx++)
00983         FOR_EACH_LEVEL(j, reverse)
00984         {
00985             if(GlobalUtil::_TruncateMethod && GlobalUtil::_FeatureCountThreshold > 0 && _featureNum > GlobalUtil::_FeatureCountThreshold) continue;
00986 
00987                         GenerateFeatureList(i, j); 
00988                         if(GlobalUtil::_timingO)        
00989             {
00990                 int idx = i * param._dog_level_num + j;
00991                 std::cout<< _levelFeatureNum[idx] <<"\t";
00992                 ocount += _levelFeatureNum[idx];
00993             }
00994                 }
00995                 if(GlobalUtil::_timingO)
00996                 {       
00997                         t2 = CLOCK(); 
00998                         std::cout << "| \t" << int(ocount) << " :\t(" << (t2 - t1) << ")\n";
00999                 }
01000         }
01001         if(GlobalUtil::_timingS)glFinish();
01002         if(GlobalUtil::_verbose)
01003         {
01004                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
01005         }
01006 }
01007 
01008 
01009 void PyramidGL::GenerateFeatureDisplayVBO()
01010 {
01011         //use a big VBO to save all the SIFT box vertices
01012         int w, h, esize; GLint bsize;
01013         int nvbo = _octave_num * param._dog_level_num;
01014         //initialize the vbos
01015         if(_featureDisplayVBO==NULL)
01016         {
01017                 _featureDisplayVBO = new GLuint[nvbo];
01018                 glGenBuffers( nvbo, _featureDisplayVBO );
01019     }
01020     if(_featurePointVBO == NULL)
01021     {
01022                 _featurePointVBO = new GLuint[nvbo];
01023                 glGenBuffers(nvbo, _featurePointVBO);
01024         }
01025 
01026         FrameBufferObject fbo;
01027         glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
01028         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
01029         glActiveTexture(GL_TEXTURE0);
01030         
01031         //
01032         GLTexImage & tempTex = *_descriptorTex;
01033         //
01034         for(int i = 0, idx = 0; i < _octave_num; i++)
01035         {
01036                 for(int j = 0; j < param._dog_level_num; j ++, idx++)
01037                 {
01038                         GLTexImage * ftex  = _featureTex + idx;
01039 
01040                         if(_levelFeatureNum[idx]<=0)continue;
01041                         //box display vbo
01042                         int count = _levelFeatureNum[idx]* 10;
01043                         GetAlignedStorageSize(count, _alignment, w, h);
01044                         w = (int)ceil(double(count)/ h);
01045 
01046                         //input
01047                         fbo.BindFBO();
01048                         ftex->BindTex();
01049 
01050                         //output
01051                         tempTex.AttachToFBO(0);
01052                         GlobalUtil::FitViewPort(w, h);
01053                         //shader
01054                         ShaderMan::UseShaderGenVBO(     (float)ftex->GetImgWidth(),  (float) w, 
01055                                 param.GetLevelSigma(j + param._level_min + 1));
01056                         GLTexImage::DrawQuad(0,  (float)w, 0, (float)h);
01057                 
01058                         //
01059                         glBindBuffer(GL_PIXEL_PACK_BUFFER_ARB, _featureDisplayVBO[ idx]);
01060                         glGetBufferParameteriv(GL_PIXEL_PACK_BUFFER_ARB, GL_BUFFER_SIZE, &bsize);
01061                         esize = w*h * sizeof(float)*4;
01062                         //increase size when necessary
01063                         if(bsize < esize) glBufferData(GL_PIXEL_PACK_BUFFER_ARB, esize*3/2,     NULL, GL_STATIC_DRAW_ARB);
01064                         
01065                         //read back if we have enough buffer    
01066                         glGetBufferParameteriv(GL_PIXEL_PACK_BUFFER_ARB, GL_BUFFER_SIZE, &bsize);
01067                         if(bsize >= esize)      glReadPixels(0, 0, w, h, GL_RGBA, GL_FLOAT, 0);
01068                         else glBufferData(GL_PIXEL_PACK_BUFFER_ARB, 0,  NULL, GL_STATIC_DRAW_ARB);
01069                         glBindBuffer(GL_PIXEL_PACK_BUFFER_ARB, 0);
01070 
01071 
01072                         //copy the texture into vbo
01073                         fbo.BindFBO();
01074                         tempTex.AttachToFBO(0);
01075 
01076                         ftex->BindTex();
01077                         ftex->FitTexViewPort();
01078                         ShaderMan::UseShaderCopyKeypoint();
01079                         ftex->DrawQuad();
01080 
01081                         glBindBuffer(GL_PIXEL_PACK_BUFFER_ARB,  _featurePointVBO[ idx]);
01082                         glGetBufferParameteriv(GL_PIXEL_PACK_BUFFER_ARB, GL_BUFFER_SIZE, &bsize);
01083                         esize = ftex->GetImgHeight() * ftex->GetImgWidth()*sizeof(float) *4;
01084 
01085                         //increase size when necessary
01086                         if(bsize < esize)       glBufferData(GL_PIXEL_PACK_BUFFER_ARB, esize*3/2 ,      NULL, GL_STATIC_DRAW_ARB);
01087 
01088                         //read back if we have enough buffer
01089                         glGetBufferParameteriv(GL_PIXEL_PACK_BUFFER_ARB, GL_BUFFER_SIZE, &bsize);
01090                         if(bsize >= esize) glReadPixels(0, 0, ftex->GetImgWidth(), ftex->GetImgHeight(), GL_RGBA, GL_FLOAT, 0);
01091                         else  glBufferData(GL_PIXEL_PACK_BUFFER_ARB, 0, NULL, GL_STATIC_DRAW_ARB);
01092 
01093                         glBindBuffer(GL_PIXEL_PACK_BUFFER_ARB, 0);
01094                         
01095                 }
01096         }
01097         glReadBuffer(GL_NONE);
01098         glFinish();
01099 
01100 }
01101 
01102 
01103 
01104 
01105 
01106 void PyramidNaive::GetFeatureOrientations()
01107 {
01108         GLTexImage * gtex;
01109         GLTexImage * stex = NULL;
01110         GLTexImage * ftex = _featureTex;
01111         GLTexImage * otex = _orientationTex;
01112         int sid = 0; 
01113         int * count      = _levelFeatureNum;
01114         float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num);
01115         FrameBufferObject fbo;
01116         if(_orientationTex)
01117         {
01118                 GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT };
01119                 glDrawBuffers(2, buffers);
01120         }else
01121         {
01122                 glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
01123         }
01124         for(int i = 0; i < _octave_num; i++)
01125         {
01126                 gtex = GetLevelTexture(i+_octave_min, param._level_min + 1);
01127                 if(GlobalUtil::_SubpixelLocalization || GlobalUtil::_KeepExtremumSign)
01128                         stex = GetBaseLevel(i+_octave_min, DATA_KEYPOINT) + 2;
01129 
01130                 for(int j = 0; j < param._dog_level_num; j++, ftex++, otex++, count++, gtex++, stex++)
01131                 {
01132                         if(*count<=0)continue;
01133 
01134                         sigma = param.GetLevelSigma(j+param._level_min+1);
01135 
01136                         //
01137                         ftex->FitTexViewPort();
01138 
01139                         glActiveTexture(GL_TEXTURE0);
01140                         ftex->BindTex();
01141                         glActiveTexture(GL_TEXTURE1);
01142                         gtex->BindTex();
01143                         //
01144                         ftex->AttachToFBO(0);
01145                         if(_orientationTex)             otex->AttachToFBO(1);
01146                         if(!_existing_keypoints && (GlobalUtil::_SubpixelLocalization|| GlobalUtil::_KeepExtremumSign))
01147                         {
01148                                 glActiveTexture(GL_TEXTURE2);
01149                                 stex->BindTex();
01150                                 sid = * stex;
01151                         }
01152                         ShaderMan::UseShaderOrientation(gtex->GetTexID(),
01153                                 gtex->GetImgWidth(), gtex->GetImgHeight(),
01154                                 sigma, sid, sigma_step, _existing_keypoints);
01155                         ftex->DrawQuad();
01156         //              glFinish();
01157                         
01158                 }
01159         }
01160 
01161         GLTexImage::UnbindMultiTex(3);
01162         if(GlobalUtil::_timingS)glFinish();
01163 
01164         if(_orientationTex)     fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
01165 
01166 }
01167 
01168 
01169 
01170 //to compare with GPU feature list generation
01171 void PyramidNaive::GenerateFeatureListCPU()
01172 {
01173 
01174         FrameBufferObject fbo;
01175         _featureNum = 0;
01176         GLTexImage * tex = GetBaseLevel(_octave_min);
01177         float * mem = new float [tex->GetTexWidth()*tex->GetTexHeight()];
01178         vector<float> list;
01179         int idx = 0;
01180         for(int i = 0; i < _octave_num; i++)
01181         {
01182                 for(int j = 0; j < param._dog_level_num; j++, idx++)
01183                 {
01184                         tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + j + 2;
01185                         tex->BindTex();
01186                         glGetTexImage(GlobalUtil::_texTarget, 0, GL_RED, GL_FLOAT, mem);
01187                         //tex->AttachToFBO(0);
01188                         //tex->FitTexViewPort();
01189                         //glReadPixels(0, 0, tex->GetTexWidth(), tex->GetTexHeight(), GL_RED, GL_FLOAT, mem);
01190                         //
01191                         //make a list of 
01192                         list.resize(0);
01193                         float * p = mem;
01194                         int fcount = 0 ;
01195                         for(int k = 0; k < tex->GetTexHeight(); k++)
01196                         {
01197                                 for( int m = 0; m < tex->GetTexWidth(); m ++, p++)
01198                                 {
01199                                         if(*p==0)continue;
01200                                         if(m ==0 || k ==0 || k >= tex->GetImgHeight() -1 || m >= tex->GetImgWidth() -1 ) continue;
01201                                         list.push_back(m+0.5f);
01202                                         list.push_back(k+0.5f);
01203                                         list.push_back(0);
01204                                         list.push_back(1);
01205                                         fcount ++;
01206 
01207 
01208                                 }
01209                         }
01210                         if(fcount==0)continue;
01211 
01212 
01213                         
01214                         GLTexImage * ftex = _featureTex+idx;
01215                         _levelFeatureNum[idx] = (fcount);
01216                         SetLevelFeatureNum(idx, fcount);
01217 
01218                         _featureNum += (fcount);
01219 
01220 
01221                         int fw = ftex->GetImgWidth();
01222                         int fh = ftex->GetImgHeight();
01223 
01224                         list.resize(4*fh*fw);
01225 
01226                         ftex->BindTex();
01227                         ftex->AttachToFBO(0);
01228         //              glTexImage2D(GlobalUtil::_texTarget, 0, GlobalUtil::_iTexFormat, fw, fh, 0, GL_BGRA, GL_FLOAT, &list[0]);
01229                         glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, fw, fh, GL_RGBA, GL_FLOAT, &list[0]);
01230                         //
01231                 }
01232         }
01233         GLTexImage::UnbindTex();
01234         delete mem;
01235         if(GlobalUtil::_verbose)
01236         {
01237                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
01238         }
01239 }
01240 
01241 #define FEATURELIST_USE_PBO
01242 
01243 void PyramidGL::ReshapeFeatureListCPU()
01244 {
01245         //make a compact feature list, each with only one orientation
01246         //download orientations and the featue list
01247         //reshape it and upload it
01248 
01249         FrameBufferObject fbo;
01250         int i, szmax =0, sz;
01251         int n = param._dog_level_num*_octave_num;
01252         for( i = 0; i < n; i++)
01253         {
01254                 sz = _featureTex[i].GetImgWidth() * _featureTex[i].GetImgHeight();
01255                 if(sz > szmax ) szmax = sz;
01256         }
01257         float * buffer = new float[szmax*24];
01258         float * buffer1 = buffer;
01259         float * buffer2 = buffer + szmax*4;
01260         float * buffer3 = buffer + szmax*8;
01261 
01262         glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
01263         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
01264 
01265 #ifdef FEATURELIST_USE_PBO
01266     GLuint ListUploadPBO;
01267     glGenBuffers(1, &ListUploadPBO); 
01268         glBindBuffer(GL_PIXEL_UNPACK_BUFFER_ARB,  ListUploadPBO);
01269         glBufferData(GL_PIXEL_UNPACK_BUFFER_ARB, szmax * 8 * sizeof(float),     NULL, GL_STREAM_DRAW);
01270 #endif
01271 
01272         _featureNum = 0;
01273 
01274 #ifdef NO_DUPLICATE_DOWNLOAD
01275         const double twopi = 2.0*3.14159265358979323846;
01276         _keypoint_buffer.resize(0);
01277         float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
01278         if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
01279         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
01280 #endif
01281 
01282         for(i = 0; i < n; i++)
01283         {
01284                 if(_levelFeatureNum[i]==0)continue;
01285 
01286                 _featureTex[i].AttachToFBO(0);
01287                 _featureTex[i].FitTexViewPort();
01288                 glReadPixels(0, 0, _featureTex[i].GetImgWidth(), _featureTex[i].GetImgHeight(),GL_RGBA, GL_FLOAT, buffer1);
01289                 
01290                 int fcount =0, ocount;
01291                 float * src = buffer1;
01292                 float * orientation  = buffer2;
01293                 float * des = buffer3;
01294                 if(GlobalUtil::_OrientationPack2 == 0)
01295                 {       
01296                         //read back orientations from another texture
01297                         _orientationTex[i].AttachToFBO(0);
01298                         glReadPixels(0, 0, _orientationTex[i].GetImgWidth(), _orientationTex[i].GetImgHeight(),GL_RGBA, GL_FLOAT, buffer2);
01299                         //make the feature list
01300                         for(int j = 0; j < _levelFeatureNum[i]; j++, src+=4, orientation+=4)
01301                         {
01302                                 if(_existing_keypoints) 
01303                                 {
01304                                         des[0] = src[0];
01305                                         des[1] = src[1];
01306                                         des[2] = orientation[0];
01307                                         des[3] = src[3];                        
01308                                         fcount++;
01309                                         des += 4;
01310                                 }else
01311                                 {
01312                                         ocount = (int)src[2];
01313                                         for(int k = 0 ; k < ocount; k++, des+=4)
01314                                         {
01315                                                 des[0] = src[0];
01316                                                 des[1] = src[1];
01317                                                 des[2] = orientation[k];
01318                                                 des[3] = src[3];                        
01319                                                 fcount++;
01320                                         }
01321                                 }
01322                         }
01323                 }else
01324                 {
01325                         _featureTex[i].DetachFBO(0);
01326                         const static double factor  = 2.0*3.14159265358979323846/65535.0;
01327                         for(int j = 0; j < _levelFeatureNum[i]; j++, src+=4)
01328                         {
01329                                 unsigned short * orientations = (unsigned short*) (&src[2]);
01330                                 if(_existing_keypoints) 
01331                                 {
01332                                         des[0] = src[0];
01333                                         des[1] = src[1];
01334                                         des[2] = float( factor* orientations[0]);
01335                                         des[3] = src[3];                        
01336                                         fcount++;
01337                                         des += 4;
01338                                 }else
01339                                 {
01340                                         if(orientations[0] != 65535)
01341                                         {
01342                                                 des[0] = src[0];
01343                                                 des[1] = src[1];
01344                                                 des[2] = float( factor* orientations[0]);
01345                                                 des[3] = src[3];                        
01346                                                 fcount++;
01347                                                 des += 4;
01348 
01349                                                 if(orientations[1] != 65535)
01350                                                 {
01351                                                         des[0] = src[0];
01352                                                         des[1] = src[1];
01353                                                         des[2] = float(factor* orientations[1]);
01354                                                         des[3] = src[3];                        
01355                                                         fcount++;
01356                                                         des += 4;
01357                                                 }
01358                                         }
01359                                 }
01360                         }
01361                 }
01362 
01363                 //texture size --------------
01364                 SetLevelFeatureNum(i, fcount);
01365                 int nfw = _featureTex[i].GetImgWidth();
01366                 int nfh = _featureTex[i].GetImgHeight();
01367                 int sz = nfh * nfw;
01368                 if(sz > fcount) memset(des, 0, sizeof(float) * (sz - fcount) * 4);
01369 
01370 #ifndef FEATURELIST_USE_PBO
01371                 _featureTex[i].BindTex();
01372                 glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, nfw, nfh, GL_RGBA, GL_FLOAT, buffer3);
01373                 _featureTex[i].UnbindTex();
01374 #else
01375         float* mem = (float*) glMapBuffer(GL_PIXEL_UNPACK_BUFFER_ARB,  GL_WRITE_ONLY);
01376         memcpy(mem, buffer3,  sz * 4 * sizeof(float) );
01377         glUnmapBuffer(GL_PIXEL_UNPACK_BUFFER_ARB);
01378                 _featureTex[i].BindTex();
01379                 glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, nfw, nfh, GL_RGBA, GL_FLOAT, 0);
01380                 _featureTex[i].UnbindTex();
01381 #endif
01382 
01383 #ifdef NO_DUPLICATE_DOWNLOAD
01384                 if(fcount > 0)
01385                 {
01386                         float oss = os * (1 << (i / param._dog_level_num));
01387                         _keypoint_buffer.resize((_featureNum + fcount) * 4);
01388                         float* ds = &_keypoint_buffer[_featureNum * 4];
01389                         float* fs = buffer3;
01390                         for(int k = 0;  k < fcount; k++, ds+=4, fs+=4)
01391                         {
01392                                 ds[0] = oss*(fs[0]-0.5f) + offset;      //x
01393                                 ds[1] = oss*(fs[1]-0.5f) + offset;      //y
01394                                 ds[3] = (float)fmod(twopi-fs[2], twopi);        //orientation, mirrored
01395                                 ds[2] = oss*fs[3];  //scale
01396                         }
01397                 }
01398 #endif
01399                 _levelFeatureNum[i] = fcount;
01400                 _featureNum += fcount;
01401         }
01402 
01403         delete[] buffer;
01404         if(GlobalUtil::_verbose)
01405         {
01406                 std::cout<<"#Features MO:\t"<<_featureNum<<endl;
01407         }
01409 #ifdef FEATURELIST_USE_PBO
01410         glBindBuffer(GL_PIXEL_UNPACK_BUFFER_ARB,  0);
01411     glDeleteBuffers(1, &ListUploadPBO);
01412 #endif
01413 }
01414 
01415 
01416 
01417 inline void PyramidGL::SetLevelFeatureNum(int idx, int fcount)
01418 {
01419         int fw, fh;
01420         GLTexImage * ftex = _featureTex + idx;
01421         //set feature texture size. normally fh will be one
01422         GetTextureStorageSize(fcount, fw, fh);
01423         if(fcount >  ftex->GetTexWidth()*ftex->GetTexHeight())
01424         {
01425         if(GlobalUtil::_timingS && GlobalUtil::_verbose)
01426                         std::cout<<"Too many features, reallocate texture\n";
01427 
01428                 ftex->InitTexture(fw, fh, 0);
01429                 if(_orientationTex)                     _orientationTex[idx].InitTexture(fw, fh, 0);
01430 
01431         }
01432         if(GlobalUtil::_NarrowFeatureTex)
01433                 fh = fcount ==0? 0:(int)ceil(double(fcount)/fw);
01434         else
01435                 fw = fcount ==0? 0:(int)ceil(double(fcount)/fh);
01436         ftex->SetImageSize(fw, fh);
01437         if(_orientationTex)             _orientationTex[idx].SetImageSize(fw, fh);
01438 }
01439 
01440 void PyramidGL::CleanUpAfterSIFT()
01441 {
01442         GLTexImage::UnbindMultiTex(3);
01443         ShaderMan::UnloadProgram();
01444         FrameBufferObject::DeleteGlobalFBO();
01445         GlobalUtil::CleanupOpenGL();
01446 }
01447 
01448 void PyramidNaive::GetSimplifiedOrientation()
01449 {
01450         //
01451         int idx = 0;
01452 //      int n = _octave_num  * param._dog_level_num;
01453         float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num); 
01454         GLTexImage * ftex = _featureTex;
01455 
01456         FrameBufferObject fbo;
01457         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
01458         for(int i = 0; i < _octave_num; i++)
01459         {
01460                 GLTexImage *gtex = GetLevelTexture(i+_octave_min, 2+param._level_min);
01461                 for(int j = 0; j < param._dog_level_num; j++, ftex++,  gtex++, idx ++)
01462                 {
01463                         if(_levelFeatureNum[idx]<=0)continue;
01464                         sigma = param.GetLevelSigma(j+param._level_min+1);
01465 
01466                         //
01467                         ftex->AttachToFBO(0);
01468                         ftex->FitTexViewPort();
01469 
01470                         glActiveTexture(GL_TEXTURE0);
01471                         ftex->BindTex();
01472                         glActiveTexture(GL_TEXTURE1);
01473                         gtex->BindTex();
01474 
01475                         ShaderMan::UseShaderSimpleOrientation(gtex->GetTexID(), sigma, sigma_step);
01476                         ftex->DrawQuad();
01477                 }
01478         }
01479 
01480         GLTexImage::UnbindMultiTex(2);
01481 
01482 }
01483 
01484 
01485 #ifdef USE_SSE_FOR_SIFTGPU
01486         static inline float dotproduct_128d(float * p)
01487         {
01488                 float z = 0.0f;
01489                 __m128 sse =_mm_load_ss(&z);
01490                 float* pf = (float*) (&sse);
01491                 for( int i = 0; i < 32; i++, p+=4)
01492                 {
01493                         __m128 ps = _mm_loadu_ps(p);
01494                         sse = _mm_add_ps(sse,  _mm_mul_ps(ps, ps));
01495                 }
01496                 return pf[0] + pf[1] + pf[2] + pf[3];
01497 
01498         }
01499         static inline void multiply_and_truncate_128d(float* p, float m)
01500         {
01501                 float z = 0.2f;
01502                 __m128 t = _mm_load_ps1(&z);
01503                 __m128 r = _mm_load_ps1(&m);
01504                 for(int i = 0; i < 32; i++, p+=4)
01505                 {
01506                         __m128 ps = _mm_loadu_ps(p);
01507                         _mm_storeu_ps(p, _mm_min_ps(_mm_mul_ps(ps, r), t));
01508                 }
01509         }
01510         static inline void multiply_128d(float* p, float m)
01511         {
01512                 __m128 r = _mm_load_ps1(&m);
01513                 for(int i = 0; i < 32; i++, p+=4)
01514                 {
01515                         __m128 ps = _mm_loadu_ps(p);
01516                         _mm_storeu_ps(p, _mm_mul_ps(ps, r));
01517                 }
01518         }
01519 #endif
01520 
01521 
01522 inline void PyramidGL::NormalizeDescriptor(int num, float*pd)
01523 {
01524 
01525 #ifdef USE_SSE_FOR_SIFTGPU
01526         for(int k = 0; k < num; k++, pd +=128)
01527         {
01528                 float sq;
01529                 //normalize and truncate to .2
01530                 sq = dotproduct_128d(pd);               sq = 1.0f / sqrtf(sq);
01531                 multiply_and_truncate_128d(pd, sq);
01532 
01533                 //renormalize
01534                 sq = dotproduct_128d(pd);               sq = 1.0f / sqrtf(sq);
01535                 multiply_128d(pd, sq);
01536         }
01537 #else
01538         //descriptor normalization runs on cpu for OpenGL implemenations
01539         for(int k = 0; k < num; k++, pd +=128)
01540         {
01541                 int v;
01542                 float* ppd, sq = 0;
01543                 //int v;
01544                 //normalize
01545                 ppd = pd;
01546                 for(v = 0 ; v < 128; v++, ppd++)        sq += (*ppd)*(*ppd);
01547                 sq = 1.0f / sqrtf(sq);
01548                 //truncate to .2
01549                 ppd = pd;
01550                 for(v = 0; v < 128; v ++, ppd++)        *ppd = min(*ppd*sq, 0.2f);
01551 
01552                 //renormalize
01553                 ppd = pd; sq = 0;
01554                 for(v = 0; v < 128; v++, ppd++) sq += (*ppd)*(*ppd);
01555                 sq = 1.0f / sqrtf(sq);
01556 
01557                 ppd = pd;
01558                 for(v = 0; v < 128; v ++, ppd++)        *ppd = *ppd*sq;
01559         }
01560 
01561 #endif
01562 }
01563 
01564 inline void PyramidGL::InterlaceDescriptorF2(int w, int h, float* buf, float* pd, int step)
01565 {
01566         /*
01567         if(GlobalUtil::_DescriptorPPR == 8)
01568         {
01569                 const int dstep = w * 128;
01570                 float* pp1 = buf;
01571                 float* pp2 = buf + step;
01572 
01573                 for(int u = 0; u < h ; u++, pd+=dstep)
01574                 {
01575                         int v; 
01576                         float* ppd = pd;
01577                         for(v= 0; v < w; v++)
01578                         {
01579                                 for(int t = 0; t < 8; t++)
01580                                 {
01581                                         *ppd++ = *pp1++;*ppd++ = *pp1++;*ppd++ = *pp1++;*ppd++ = *pp1++;
01582                                         *ppd++ = *pp2++;*ppd++ = *pp2++;*ppd++ = *pp2++;*ppd++ = *pp2++;
01583                                 }
01584                                 ppd += 64;
01585                         }
01586                         ppd = pd + 64;
01587                         for(v= 0; v < w; v++)
01588                         {
01589                                 for(int t = 0; t < 8; t++)
01590                                 {
01591                                         *ppd++ = *pp1++;*ppd++ = *pp1++;*ppd++ = *pp1++;*ppd++ = *pp1++;
01592                                         *ppd++ = *pp2++;*ppd++ = *pp2++;*ppd++ = *pp2++;*ppd++ = *pp2++;
01593                                 }
01594                                 ppd += 64;
01595                         }
01596                 }
01597 
01598         }else */
01599         if(GlobalUtil::_DescriptorPPR == 8)
01600         {
01601                 //interlace
01602                 for(int k = 0; k < 2; k++)
01603                 {
01604                         float* pp = buf + k * step;
01605                         float* ppd = pd + k * 4;
01606                         for(int u = 0; u < h ; u++)
01607                         {
01608                                 int v; 
01609                                 for(v= 0; v < w; v++)
01610                                 {
01611                                         for(int t = 0; t < 8; t++)
01612                                         {
01613                                                 ppd[0] = pp[0];
01614                                                 ppd[1] = pp[1];
01615                                                 ppd[2] = pp[2];
01616                                                 ppd[3] = pp[3];
01617                                                 ppd += 8;
01618                                                 pp+= 4;
01619                                         }
01620                                         ppd += 64;
01621                                 }
01622                                 ppd += ( 64 - 128 * w );
01623                                 for(v= 0; v < w; v++)
01624                                 {
01625                                         for(int t = 0; t < 8; t++)
01626                                         {
01627                                                 ppd[0] = pp[0];
01628                                                 ppd[1] = pp[1];
01629                                                 ppd[2] = pp[2];
01630                                                 ppd[3] = pp[3];
01631 
01632                                                 ppd += 8;
01633                                                 pp+= 4;
01634                                         }
01635                                         ppd += 64;
01636                                 }
01637                                 ppd -=64;
01638                         }
01639                 }
01640         }else if(GlobalUtil::_DescriptorPPR == 4)
01641         {
01642 
01643         }
01644 
01645 
01646 
01647 }
01648 void PyramidGL::GetFeatureDescriptors()
01649 {
01650         //descriptors...
01651         float sigma;
01652         int idx, i, j, k, w, h;
01653         int ndf = 32 / GlobalUtil::_DescriptorPPT; //number of textures
01654         int block_width = GlobalUtil::_DescriptorPPR;
01655         int block_height = GlobalUtil::_DescriptorPPT/GlobalUtil::_DescriptorPPR;
01656         float* pd =  &_descriptor_buffer[0], * pbuf  = NULL;
01657         vector<float>read_buffer, descriptor_buffer2;
01658 
01659         //use another buffer, if we need to re-order the descriptors
01660         if(_keypoint_index.size() > 0)
01661         {
01662                 descriptor_buffer2.resize(_descriptor_buffer.size());
01663                 pd = &descriptor_buffer2[0];
01664         }
01665         FrameBufferObject fbo;
01666 
01667         GLTexImage * gtex, *otex, * ftex;
01668         GLenum buffers[8] = { 
01669                 GL_COLOR_ATTACHMENT0_EXT,               GL_COLOR_ATTACHMENT1_EXT ,
01670                 GL_COLOR_ATTACHMENT2_EXT,               GL_COLOR_ATTACHMENT3_EXT ,
01671                 GL_COLOR_ATTACHMENT4_EXT,               GL_COLOR_ATTACHMENT5_EXT ,
01672                 GL_COLOR_ATTACHMENT6_EXT,               GL_COLOR_ATTACHMENT7_EXT ,
01673         };
01674 
01675         glDrawBuffers(ndf, buffers);
01676         glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
01677 
01678 
01679         for( i = 0, idx = 0, ftex = _featureTex; i < _octave_num; i++)
01680         {
01681                 gtex = GetBaseLevel(i + _octave_min, DATA_GRAD) + 1;
01682                 otex = GetBaseLevel(i + _octave_min, DATA_ROT)  + 1;
01683                 for( j = 0; j < param._dog_level_num; j++, ftex++, idx++, gtex++, otex++)
01684                 {
01685                         if(_levelFeatureNum[idx]==0)continue;
01686 
01687             sigma = IsUsingRectDescription()?  0 : param.GetLevelSigma(j+param._level_min+1);
01688                         int count = _levelFeatureNum[idx] * block_width;
01689                         GetAlignedStorageSize(count, block_width, w, h);
01690                         h = ((int)ceil(double(count) / w)) * block_height;
01691 
01692                         //not enought space for holding the descriptor data
01693                         if(w > _descriptorTex[0].GetTexWidth() || h > _descriptorTex[0].GetTexHeight())
01694                         {
01695                                 for(k = 0; k < ndf; k++)_descriptorTex[k].InitTexture(w, h);
01696                         }
01697                         for(k = 0; k < ndf; k++)        _descriptorTex[k].AttachToFBO(k);
01698                         GlobalUtil::FitViewPort(w, h);
01699                         glActiveTexture(GL_TEXTURE0);
01700                         ftex->BindTex();
01701                         glActiveTexture(GL_TEXTURE1);
01702                         gtex->BindTex();
01703                         if(otex!=gtex)
01704                         {
01705                                 glActiveTexture(GL_TEXTURE2);
01706                                 otex->BindTex();
01707                         }
01708 
01709                         ShaderMan::UseShaderDescriptor(gtex->GetTexID(), otex->GetTexID(), 
01710                                 w, ftex->GetImgWidth(), gtex->GetImgWidth(), gtex->GetImgHeight(), sigma);
01711                         GLTexImage::DrawQuad(0, (float)w, 0, (float)h);
01712 
01713                          //read back float format descriptors and do normalization on CPU
01714                         int step = w*h*4;
01715                         if((unsigned int)step*ndf > read_buffer.size())
01716                         {
01717                                 read_buffer.resize(ndf*step);
01718                         }
01719                         pbuf = &read_buffer[0];
01720                         
01721                         //read back
01722                         for(k = 0; k < ndf; k++, pbuf+=step)
01723                         {
01724                                 glReadBuffer(GL_COLOR_ATTACHMENT0_EXT + k);
01725                                 glReadPixels(0, 0, w, h, GL_RGBA, GL_FLOAT, pbuf);
01726                         }
01727         
01728                         //the following two steps run on cpu, so better cpu better speed
01729                         //and release version can be a lot faster than debug version
01730                         //interlace data on the two texture to get the descriptor
01731                         InterlaceDescriptorF2(w / block_width, h / block_height, &read_buffer[0], pd, step);
01732                         
01733                         //need to do normalization
01734                         //the new version uses SSE to speed up this part
01735                         if(GlobalUtil::_NormalizedSIFT) NormalizeDescriptor(_levelFeatureNum[idx], pd);
01736 
01737                         pd += 128*_levelFeatureNum[idx];
01738                         glReadBuffer(GL_NONE);
01739                 }
01740         }
01741 
01742 
01743         //finally, put the descriptor back to their original order for existing keypoint list.
01744         if(_keypoint_index.size() > 0)
01745         {
01746                 for(i = 0; i < _featureNum; ++i)
01747                 {
01748                         int index = _keypoint_index[i];
01749                         memcpy(&_descriptor_buffer[index*128], &descriptor_buffer2[i*128], 128 * sizeof(float));
01750                 }
01751         }
01752 
01754         GLTexImage::UnbindMultiTex(3); 
01755         glDrawBuffer(GL_NONE);
01756         ShaderMan::UnloadProgram();
01757         if(GlobalUtil::_timingS)glFinish();
01758         for(i = 0; i < ndf; i++) fbo.UnattachTex(GL_COLOR_ATTACHMENT0_EXT +i);
01759 
01760 }
01761 
01762 
01763 void PyramidGL::DownloadKeypoints()
01764 {
01765         const double twopi = 2.0*3.14159265358979323846;
01766         int idx = 0;
01767         float * buffer = &_keypoint_buffer[0];
01768         vector<float> keypoint_buffer2;
01769         //use a different keypoint buffer when processing with an exisint features list
01770         //without orientation information. 
01771         if(_keypoint_index.size() > 0)
01772         {
01773                 keypoint_buffer2.resize(_keypoint_buffer.size());
01774                 buffer = &keypoint_buffer2[0];
01775         }
01776         float * p = buffer, *ps, sigma;
01777         GLTexImage * ftex = _featureTex;
01778         FrameBufferObject fbo;
01779         ftex->FitRealTexViewPort();
01781         float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
01782         if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
01783         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
01785         for(int i = 0; i < _octave_num; i++, os *= 2.0f)
01786         {
01787                 
01788                 for(int j = 0; j  < param._dog_level_num; j++, idx++, ftex++)
01789                 {
01790 
01791                         if(_levelFeatureNum[idx]>0)
01792                         {       
01793                                 ftex->AttachToFBO(0);
01794                                 glReadPixels(0, 0, ftex->GetImgWidth(), ftex->GetImgHeight(),GL_RGBA, GL_FLOAT, p);
01795                                 ps = p;
01796                                 for(int k = 0;  k < _levelFeatureNum[idx]; k++, ps+=4)
01797                                 {
01798                                         ps[0] = os*(ps[0]-0.5f) + offset;       //x
01799                                         ps[1] = os*(ps[1]-0.5f) + offset;       //y
01800                                         sigma = os*ps[3]; 
01801                                         ps[3] = (float)fmod(twopi-ps[2], twopi);        //orientation, mirrored
01802                                         ps[2] = sigma;  //scale
01803                                 }
01804                                 p+= 4* _levelFeatureNum[idx];
01805                         }
01806                 }
01807         }
01808 
01809         //put the feature into their original order
01810 
01811         if(_keypoint_index.size() > 0)
01812         {
01813                 for(int i = 0; i < _featureNum; ++i)
01814                 {
01815                         int index = _keypoint_index[i];
01816                         memcpy(&_keypoint_buffer[index*4], &keypoint_buffer2[i*4], 4 * sizeof(float));
01817                 }
01818         }
01819 }
01820 
01821 
01822 void PyramidGL::GenerateFeatureListTex()
01823 {
01824         //generate feature list texture from existing keypoints
01825         //do feature sorting in the same time?
01826 
01827         FrameBufferObject fbo;
01828         vector<float> list;
01829         int idx = 0;
01830         const double twopi = 2.0*3.14159265358979323846;
01831         float sigma_half_step = powf(2.0f, 0.5f / param._dog_level_num);
01832         float octave_sigma = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
01833         float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f; 
01834         if(_down_sample_factor>0) octave_sigma *= float(1<<_down_sample_factor); 
01835 
01836     
01837     std::fill(_levelFeatureNum, _levelFeatureNum + _octave_num * param._dog_level_num, 0);
01838 
01839         _keypoint_index.resize(0); // should already be 0
01840         for(int i = 0; i < _octave_num; i++, octave_sigma*= 2.0f)
01841         {
01842                 for(int j = 0; j < param._dog_level_num; j++, idx++)
01843                 {
01844                         list.resize(0);
01845                         float level_sigma = param.GetLevelSigma(j + param._level_min + 1) * octave_sigma;
01846                         float sigma_min = level_sigma / sigma_half_step;
01847                         float sigma_max = level_sigma * sigma_half_step;
01848                         int fcount = 0 ;
01849                         for(int k = 0; k < _featureNum; k++)
01850                         {
01851                                 float * key = &_keypoint_buffer[k*4];
01852                 float sigmak = key[2]; 
01853 
01855                 if(IsUsingRectDescription()) sigmak = min(key[2], key[3]) / 12.0f; 
01856 
01857                 if( (sigmak >= sigma_min && sigmak < sigma_max)
01858                                         ||(sigmak < sigma_min && i ==0 && j == 0)
01859                                         ||(sigmak > sigma_max && j == param._dog_level_num - 1&& 
01860                             (i == _octave_num -1  || GlobalUtil::_KeyPointListForceLevel0)))
01861                                 {
01862                                         //add this keypoint to the list
01863                                         list.push_back((key[0] - offset) / octave_sigma + 0.5f);
01864                                         list.push_back((key[1] - offset) / octave_sigma + 0.5f);
01865                     if(IsUsingRectDescription())
01866                     {
01867                         list.push_back(key[2] / octave_sigma);
01868                         list.push_back(key[3] / octave_sigma);
01869                     }else
01870                     {
01871                                             list.push_back((float)fmod(twopi-key[3], twopi));
01872                                             list.push_back(key[2] / octave_sigma);
01873                     }
01874                                         fcount ++;
01875                                         //save the index of keypoints
01876                                         _keypoint_index.push_back(k);
01877                                 }
01878                         }
01879 
01880                         _levelFeatureNum[idx] = fcount;
01881                         if(fcount==0)continue;
01882                         GLTexImage * ftex = _featureTex+idx;
01883 
01884                         SetLevelFeatureNum(idx, fcount);
01885 
01886                         int fw = ftex->GetImgWidth();
01887                         int fh = ftex->GetImgHeight();
01888 
01889                         list.resize(4*fh*fw);
01890 
01891                         ftex->BindTex();
01892                         ftex->AttachToFBO(0);
01893                         glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, fw, fh, GL_RGBA, GL_FLOAT, &list[0]);
01894 
01895             if( fcount == _featureNum) _keypoint_index.resize(0);
01896                 }
01897         if( GlobalUtil::_KeyPointListForceLevel0 ) break;
01898         }
01899         GLTexImage::UnbindTex();
01900         if(GlobalUtil::_verbose)
01901         {
01902                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
01903         }
01904 }
01905 
01906 
01907 
01908 PyramidPacked::PyramidPacked(SiftParam& sp): PyramidGL(sp)
01909 {
01910         _allPyramid = NULL;
01911 }
01912 
01913 PyramidPacked::~PyramidPacked()
01914 {
01915         DestroyPyramidData();
01916 }
01917 
01918 
01919 //build the gaussian pyrmaid
01920 
01921 void PyramidPacked::BuildPyramid(GLTexInput * input)
01922 {
01923         //
01924         USE_TIMING();
01925         GLTexImage * tex, *tmp;
01926         FilterProgram ** filter;
01927         FrameBufferObject fbo;
01928 
01929         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
01930         input->FitTexViewPort();
01931 
01932         for (int i = _octave_min; i < _octave_min + _octave_num; i++)
01933         {
01934                 tex = GetBaseLevel(i);
01935                 tmp = GetBaseLevel(i, DATA_DOG) + 2; //use this as a temperory texture
01936 
01937                 
01938                 filter = ShaderMan::s_bag->f_gaussian_step;
01939 
01940                 OCTAVE_START();
01941 
01942                 if( i == _octave_min )
01943                 {
01944                         if(i < 0)       TextureUpSample(tex, input, 1<<(-i-1));                 
01945                         else        TextureDownSample(tex, input, 1<<(i+1));
01946             ShaderMan::FilterInitialImage(tex, tmp); 
01947                 }else
01948                 {
01949                         TextureDownSample(tex, GetLevelTexture(i-1, param._level_ds)); 
01950                         ShaderMan::FilterSampledImage(tex, tmp); 
01951                 }
01952                 LEVEL_FINISH();         
01953 
01954                 for(int j = param._level_min + 1; j <=  param._level_max ; j++, tex++, filter++)
01955                 {
01956                         // filtering
01957             ShaderMan::FilterImage(*filter, tex+1, tex, tmp);
01958                         LEVEL_FINISH();
01959                 }
01960 
01961                 OCTAVE_FINISH();
01962 
01963         }
01964         if(GlobalUtil::_timingS)        glFinish();
01965         UnloadProgram();        
01966 }
01967 
01968 void PyramidPacked::ComputeGradient()
01969 {
01970         
01971         //first pass, compute dog, gradient, orientation
01972         GLenum buffers[4] = { 
01973                 GL_COLOR_ATTACHMENT0_EXT,               GL_COLOR_ATTACHMENT1_EXT ,
01974                 GL_COLOR_ATTACHMENT2_EXT,               GL_COLOR_ATTACHMENT3_EXT
01975         };
01976 
01977         int i, j;
01978         double ts, t1;
01979         FrameBufferObject fbo;
01980 
01981         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
01982 
01983         for(i = _octave_min; i < _octave_min + _octave_num; i++)
01984         {
01985                 GLTexImage * gus = GetBaseLevel(i) +  1;
01986                 GLTexImage * dog = GetBaseLevel(i, DATA_DOG) +  1;
01987                 GLTexImage * grd = GetBaseLevel(i, DATA_GRAD) +  1;
01988                 GLTexImage * rot = GetBaseLevel(i, DATA_ROT) +  1;
01989                 glDrawBuffers(3, buffers);
01990                 gus->FitTexViewPort();
01991                 //compute the gradient
01992                 for(j = 0; j <  param._dog_level_num ; j++, gus++, dog++, grd++, rot++)
01993                 {
01994                         //gradient, dog, orientation
01995                         glActiveTexture(GL_TEXTURE0);
01996                         gus->BindTex();
01997                         glActiveTexture(GL_TEXTURE1);
01998                         (gus-1)->BindTex();
01999                         //output
02000                         dog->AttachToFBO(0);
02001                         grd->AttachToFBO(1);
02002                         rot->AttachToFBO(2);
02003                         ShaderMan::UseShaderGradientPass((gus-1)->GetTexID());
02004                         //compute
02005                         dog->DrawQuadMT4();
02006                 }
02007         }
02008         if(GlobalUtil::_timingS)
02009         {
02010                 glFinish();
02011                 if(GlobalUtil::_verbose)
02012                 {
02013                         t1 = CLOCK();
02014                         std::cout       <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n";
02015                 }
02016         }
02017         GLTexImage::DetachFBO(1);
02018         GLTexImage::DetachFBO(2);
02019         UnloadProgram();
02020         GLTexImage::UnbindMultiTex(3);
02021         fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
02022 }
02023 
02024 void PyramidPacked::DetectKeypointsEX()
02025 {
02026 
02027         //first pass, compute dog, gradient, orientation
02028         GLenum buffers[4] = { 
02029                 GL_COLOR_ATTACHMENT0_EXT,               GL_COLOR_ATTACHMENT1_EXT ,
02030                 GL_COLOR_ATTACHMENT2_EXT,               GL_COLOR_ATTACHMENT3_EXT
02031         };
02032 
02033         int i, j;
02034         double t0, t, ts, t1, t2;
02035         FrameBufferObject fbo;
02036 
02037         if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
02038 
02039         for(i = _octave_min; i < _octave_min + _octave_num; i++)
02040         {
02041                 GLTexImage * gus = GetBaseLevel(i) + 1;
02042                 GLTexImage * dog = GetBaseLevel(i, DATA_DOG) + 1;
02043                 GLTexImage * grd = GetBaseLevel(i, DATA_GRAD) + 1;
02044                 GLTexImage * rot = GetBaseLevel(i, DATA_ROT) + 1;
02045                 glDrawBuffers(3, buffers);
02046                 gus->FitTexViewPort();
02047                 //compute the gradient
02048                 for(j = param._level_min +1; j <=  param._level_max ; j++, gus++, dog++, grd++, rot++)
02049                 {
02050                         //gradient, dog, orientation
02051                         glActiveTexture(GL_TEXTURE0);
02052                         gus->BindTex();
02053                         glActiveTexture(GL_TEXTURE1);
02054                         (gus-1)->BindTex();
02055                         //output
02056                         dog->AttachToFBO(0);
02057                         grd->AttachToFBO(1);
02058                         rot->AttachToFBO(2);
02059                         ShaderMan::UseShaderGradientPass((gus-1)->GetTexID());
02060                         //compute
02061                         dog->DrawQuadMT4();
02062                 }
02063         }
02064         if(GlobalUtil::_timingS && GlobalUtil::_verbose)
02065         {
02066                 glFinish();
02067                 t1 = CLOCK();
02068         }
02069         GLTexImage::DetachFBO(1);
02070         GLTexImage::DetachFBO(2);
02071         //glDrawBuffers(1, buffers);
02072         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
02073         
02074 
02075         GlobalUtil::CheckErrorsGL();
02076 
02077         for ( i = _octave_min; i < _octave_min + _octave_num; i++)
02078         {
02079                 if(GlobalUtil::_timingO)
02080                 {
02081                         t0 = CLOCK();
02082                         std::cout<<"#"<<(i + _down_sample_factor)<<"\t";
02083                 }
02084                 GLTexImage * dog = GetBaseLevel(i, DATA_DOG) + 2;
02085                 GLTexImage * key = GetBaseLevel(i, DATA_KEYPOINT) +2;
02086                 key->FitTexViewPort();
02087 
02088                 for( j = param._level_min +2; j <  param._level_max ; j++, dog++, key++)
02089                 {
02090                         if(GlobalUtil::_timingL)t = CLOCK();            
02091                         key->AttachToFBO(0);
02092                         glActiveTexture(GL_TEXTURE0);
02093                         dog->BindTex();
02094                         glActiveTexture(GL_TEXTURE1);
02095                         (dog+1)->BindTex();
02096                         glActiveTexture(GL_TEXTURE2);
02097                         (dog-1)->BindTex();
02098                         ShaderMan::UseShaderKeypoint((dog+1)->GetTexID(), (dog-1)->GetTexID());
02099                         key->DrawQuadMT8();
02100                         if(GlobalUtil::_timingL)
02101                         {
02102                                 glFinish();
02103                                 std::cout<<(CLOCK()-t)<<"\t";
02104                         }
02105                 }
02106                 if(GlobalUtil::_timingO)
02107                 {
02108                         glFinish();
02109                         std::cout<<"|\t"<<(CLOCK()-t0)<<"\n";
02110                 }
02111         }
02112 
02113         if(GlobalUtil::_timingS)
02114         {
02115                 glFinish();
02116                 if(GlobalUtil::_verbose) 
02117                 {       
02118                         t2 = CLOCK();
02119                         std::cout       <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n"
02120                                                 <<"<Get Keypoints  >\t"<<(t2-t1)<<"\n";
02121                 }
02122                                                 
02123         }
02124         UnloadProgram();
02125         GLTexImage::UnbindMultiTex(3);
02126         fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
02127 }
02128 
02129 
02130 void PyramidPacked::GenerateFeatureList(int i, int j)
02131 {
02132         float fcount = 0.0f; 
02133         int hist_skip_gpu = GlobalUtil::_ListGenSkipGPU;
02134     int idx = i * param._dog_level_num + j; 
02135     int hist_level_num = _hpLevelNum - _pyramid_octave_first;
02136         GLTexImage * htex, * ftex, * tex;
02137         htex = _histoPyramidTex + hist_level_num - 1 - i;
02138         ftex = _featureTex + idx;
02139         tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2 + j;
02140 
02141 
02142         //fill zero to an extra row/col if the height/width is odd
02143         glActiveTexture(GL_TEXTURE0);
02144         tex->BindTex();
02145         htex->AttachToFBO(0);
02146         int tight = (htex->GetImgWidth() * 4 == tex->GetImgWidth() -1 && htex->GetImgHeight() *4 == tex->GetImgHeight()-1 );
02147         ShaderMan::UseShaderGenListInit(tex->GetImgWidth(), tex->GetImgHeight(), tight);
02148         htex->FitTexViewPort();
02149         //this uses the fact that no feature is on the edge.
02150         htex->DrawQuadReduction();
02151         //reduction..
02152         htex--;
02153         
02154         //this part might have problems on several GPUS
02155         //because the output of one pass is the input of the next pass
02156         //may require glFinish to make it right, but too much glFinish makes it slow
02157         for(int k = 0; k <hist_level_num - i-1 - hist_skip_gpu; k++, htex--)
02158         {
02159                 htex->AttachToFBO(0);
02160                 htex->FitTexViewPort();
02161                 (htex+1)->BindTex();
02162                 ShaderMan::UseShaderGenListHisto();
02163                 htex->DrawQuadReduction();              
02164         }
02165 
02166         if(hist_skip_gpu == 0)
02167         {               
02168                 //read back one pixel
02169                 float fn[4];
02170                 glReadPixels(0, 0, 1, 1, GL_RGBA , GL_FLOAT, fn);
02171                 fcount = (fn[0] + fn[1] + fn[2] + fn[3]);
02172                 if(fcount < 1) fcount = 0;
02173 
02174                 _levelFeatureNum[ idx] = (int)(fcount);
02175                 SetLevelFeatureNum(idx, (int)fcount);
02176 
02177                 //save  number of features
02178                 _featureNum += int(fcount);
02179 
02180                 //
02181                 if(fcount < 1.0)                                return;;
02182 
02183 
02185                 htex=  _histoPyramidTex;
02186 
02187                 htex->BindTex();
02188 
02189                 //first pass
02190                 ftex->AttachToFBO(0);
02191                 if(GlobalUtil::_MaxOrientation>1)
02192                 {
02193                         //this is very important...
02194                         ftex->FitRealTexViewPort();
02195                         glClear(GL_COLOR_BUFFER_BIT);
02196                         glFinish();
02197                 }else
02198                 {       
02199                         ftex->FitTexViewPort();
02200                         //glFinish();
02201                 }
02202 
02203 
02204                 ShaderMan::UseShaderGenListStart((float)ftex->GetImgWidth(), htex->GetTexID());
02205 
02206                 ftex->DrawQuad();
02207                 //make sure it finishes before the next step
02208                 ftex->DetachFBO(0);
02209                 //pass on each pyramid level
02210                 htex++;
02211         }else
02212         {
02213 
02214                 int tw = htex[1].GetDrawWidth(), th = htex[1].GetDrawHeight();
02215                 int fc = 0;
02216                 glReadPixels(0, 0, tw, th, GL_RGBA , GL_FLOAT, _histo_buffer);  
02217                 _keypoint_buffer.resize(0);
02218                 for(int y = 0, pos = 0; y < th; y++)
02219                 {
02220                         for(int x= 0; x < tw; x++)
02221                         {
02222                                 for(int c = 0; c < 4; c++, pos++)
02223                                 {
02224                                         int ss =  (int) _histo_buffer[pos]; 
02225                                         if(ss == 0) continue;
02226                                         float ft[4] = {2 * x + (c%2? 1.5f:  0.5f), 2 * y + (c>=2? 1.5f: 0.5f), 0, 1 };
02227                                         for(int t = 0; t < ss; t++)
02228                                         {
02229                                                 ft[2] = (float) t; 
02230                                                 _keypoint_buffer.insert(_keypoint_buffer.end(), ft, ft+4);
02231                                         }
02232                                         fc += (int)ss; 
02233                                 }
02234                         }
02235                 }
02236                 _levelFeatureNum[ idx] = fc;
02237                 SetLevelFeatureNum(idx, fc);
02238                 if(fc == 0)  return; 
02239 
02240                 fcount = (float) fc;    
02241                 _featureNum += fc;
02243                 ftex->AttachToFBO(0);
02244                 if(GlobalUtil::_MaxOrientation>1)
02245                 {
02246                         ftex->FitRealTexViewPort();
02247                         glClear(GL_COLOR_BUFFER_BIT);
02248                 }else
02249                 {                                       
02250                         ftex->FitTexViewPort();
02251             glFlush();
02252                 }
02253                 _keypoint_buffer.resize(ftex->GetDrawWidth() * ftex->GetDrawHeight()*4, 0);
02255                 glActiveTexture(GL_TEXTURE0);
02256                 ftex->BindTex();
02257                 glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, ftex->GetDrawWidth(),
02258                         ftex->GetDrawHeight(), GL_RGBA, GL_FLOAT, &_keypoint_buffer[0]);
02259                 htex += 2; 
02260         }
02261 
02262         for(int lev = 1 + hist_skip_gpu; lev < hist_level_num  - i; lev++, htex++)
02263         {
02264                 glActiveTexture(GL_TEXTURE0);
02265                 ftex->BindTex();
02266                 ftex->AttachToFBO(0);
02267                 glActiveTexture(GL_TEXTURE1);
02268                 htex->BindTex();
02269                 ShaderMan::UseShaderGenListStep(ftex->GetTexID(), htex->GetTexID());
02270                 ftex->DrawQuad();
02271                 ftex->DetachFBO(0);     
02272         }
02273 
02274         ftex->AttachToFBO(0);
02275         glActiveTexture(GL_TEXTURE1);
02276         tex->BindTex();
02277         ShaderMan::UseShaderGenListEnd(tex->GetTexID());
02278         ftex->DrawQuad();
02279         GLTexImage::UnbindMultiTex(2);
02280 
02281 }
02282 
02283 void PyramidPacked::GenerateFeatureList()
02284 {
02285         //generate the histogram pyramid
02286         FrameBufferObject fbo;
02287         glReadBuffer(GL_COLOR_ATTACHMENT0_EXT);
02288         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
02289         double t1, t2; 
02290         int ocount= 0, reverse = (GlobalUtil::_TruncateMethod == 1);
02291         _featureNum = 0;
02292 
02293         FitHistogramPyramid();
02294         //for(int i = 0, idx = 0; i < _octave_num; i++)
02295     FOR_EACH_OCTAVE(i, reverse)
02296         {
02297                 if(GlobalUtil::_timingO)
02298                 {
02299             t1= CLOCK();
02300                         ocount = 0;
02301                         std::cout<<"#"<<i+_octave_min + _down_sample_factor<<":\t";
02302                 }
02303                 //for(int j = 0; j < param._dog_level_num; j++, idx++)
02304         FOR_EACH_LEVEL(j, reverse)
02305                 {
02306             if(GlobalUtil::_TruncateMethod && GlobalUtil::_FeatureCountThreshold > 0 && _featureNum > GlobalUtil::_FeatureCountThreshold) continue;
02307             
02308             GenerateFeatureList(i, j); 
02309 
02310                         if(GlobalUtil::_timingO)
02311                         {
02312                 int idx = i * param._dog_level_num + j;
02313                 ocount += _levelFeatureNum[idx];
02314                                 std::cout<< _levelFeatureNum[idx] <<"\t";
02315                         }
02316                 }
02317                 if(GlobalUtil::_timingO)
02318                 {       
02319             t2 = CLOCK();
02320                         std::cout << "| \t" << int(ocount) << " :\t(" << (t2 - t1) << ")\n";
02321                 }
02322         }
02323         if(GlobalUtil::_timingS)glFinish();
02324         if(GlobalUtil::_verbose)
02325         {
02326                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
02327         }
02328 
02329 }
02330 
02331 void PyramidPacked::GenerateFeatureListCPU()
02332 {
02333         FrameBufferObject fbo;
02334         _featureNum = 0;
02335         GLTexImage * tex = GetBaseLevel(_octave_min);
02336         float * mem = new float [tex->GetTexWidth()*tex->GetTexHeight()*4];
02337         vector<float> list;
02338         int idx = 0;
02339         for(int i = 0; i < _octave_num; i++)
02340         {
02341                 for(int j = 0; j < param._dog_level_num; j++, idx++)
02342                 {
02343                         tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + j + 2;
02344                         tex->BindTex();
02345                         glGetTexImage(GlobalUtil::_texTarget, 0, GL_RGBA, GL_FLOAT, mem);
02346                         //tex->AttachToFBO(0);
02347                         //tex->FitTexViewPort();
02348                         //glReadPixels(0, 0, tex->GetTexWidth(), tex->GetTexHeight(), GL_RED, GL_FLOAT, mem);
02349                         //
02350                         //make a list of 
02351                         list.resize(0);
02352                         float *pl = mem;
02353                         int fcount = 0 ;
02354                         for(int k = 0; k < tex->GetDrawHeight(); k++)
02355                         {
02356                                 float * p = pl; 
02357                                 pl += tex->GetTexWidth() * 4;
02358                                 for( int m = 0; m < tex->GetDrawWidth(); m ++, p+=4)
02359                                 {
02360                                 //      if(m ==0 || k ==0 || k == tex->GetDrawHeight() -1 || m == tex->GetDrawWidth() -1) continue;
02361                                 //      if(*p == 0) continue;
02362                                         int t = ((int) fabs(p[0])) - 1;
02363                                         if(t < 0) continue;
02364                                         int xx = m + m + ( (t %2)? 1 : 0);
02365                                         int yy = k + k + ( (t <2)? 0 : 1);
02366                                         if(xx ==0 || yy == 0) continue;
02367                                         if(xx >= tex->GetImgWidth() - 1 || yy >= tex->GetImgHeight() - 1)continue;
02368                                         list.push_back(xx + 0.5f + p[1]);
02369                                         list.push_back(yy + 0.5f + p[2]);
02370                                         list.push_back(GlobalUtil::_KeepExtremumSign && p[0] < 0 ? -1.0f : 1.0f);
02371                                         list.push_back(p[3]);
02372                                         fcount ++;
02373                                 }
02374                         }
02375                         if(fcount==0)continue;
02376 
02377                         if(GlobalUtil::_timingL) std::cout<<fcount<<".";
02378                         
02379                         GLTexImage * ftex = _featureTex+idx;
02380                         _levelFeatureNum[idx] = (fcount);
02381                         SetLevelFeatureNum(idx, fcount);
02382 
02383                         _featureNum += (fcount);
02384 
02385 
02386                         int fw = ftex->GetImgWidth();
02387                         int fh = ftex->GetImgHeight();
02388 
02389                         list.resize(4*fh*fw);
02390 
02391                         ftex->BindTex();
02392                         ftex->AttachToFBO(0);
02393         //              glTexImage2D(GlobalUtil::_texTarget, 0, GlobalUtil::_iTexFormat, fw, fh, 0, GL_BGRA, GL_FLOAT, &list[0]);
02394                         glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, fw, fh, GL_RGBA, GL_FLOAT, &list[0]);
02395                         //
02396                 }
02397         }
02398         GLTexImage::UnbindTex();
02399         delete mem;
02400         if(GlobalUtil::_verbose)
02401         {
02402                 std::cout<<"#Features:\t"<<_featureNum<<"\n";
02403         }
02404 }
02405 
02406 
02407 
02408 void PyramidPacked::GetFeatureOrientations()
02409 {
02410         GLTexImage * gtex, * otex;
02411         GLTexImage * ftex = _featureTex;
02412         GLTexImage * fotex = _orientationTex; 
02413         int * count      = _levelFeatureNum;
02414         float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num);
02415 
02416 
02417         FrameBufferObject fbo;
02418         if(_orientationTex)
02419         {
02420                 GLenum buffers[] = { GL_COLOR_ATTACHMENT0_EXT, GL_COLOR_ATTACHMENT1_EXT };
02421                 glDrawBuffers(2, buffers);
02422         }else
02423         {
02424                 glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
02425         }
02426 
02427         for(int i = 0; i < _octave_num; i++)
02428         {
02429                 gtex = GetBaseLevel(i+_octave_min, DATA_GRAD) + 1;
02430                 otex = GetBaseLevel(i+_octave_min, DATA_ROT) + 1;
02431 
02432                 for(int j = 0; j < param._dog_level_num; j++, ftex++, otex++, count++, gtex++, fotex++)
02433                 {
02434                         if(*count<=0)continue;
02435 
02436                         sigma = param.GetLevelSigma(j+param._level_min+1);
02437 
02438 
02439                         ftex->FitTexViewPort();
02440 
02441                         glActiveTexture(GL_TEXTURE0);
02442                         ftex->BindTex();
02443                         glActiveTexture(GL_TEXTURE1);
02444                         gtex->BindTex();
02445                         glActiveTexture(GL_TEXTURE2);
02446                         otex->BindTex();
02447                         //
02448                         ftex->AttachToFBO(0);
02449                         if(_orientationTex)             fotex->AttachToFBO(1);
02450 
02451                         ShaderMan::UseShaderOrientation(gtex->GetTexID(),
02452                                 gtex->GetImgWidth(), gtex->GetImgHeight(),
02453                                 sigma, otex->GetTexID(), sigma_step, _existing_keypoints);
02454                         ftex->DrawQuad();
02455                 }
02456         }
02457 
02458         GLTexImage::UnbindMultiTex(3);
02459         if(GlobalUtil::_timingS)glFinish();
02460         if(_orientationTex)     fbo.UnattachTex(GL_COLOR_ATTACHMENT1_EXT);
02461 
02462 }
02463 
02464 
02465 void PyramidPacked::GetSimplifiedOrientation()
02466 {
02467         //
02468         int idx = 0;
02469 //      int n = _octave_num  * param._dog_level_num;
02470         float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num); 
02471         GLTexImage * ftex = _featureTex;
02472 
02473         FrameBufferObject fbo;
02474         glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
02475         for(int i = 0; i < _octave_num; i++)
02476         {
02477                 GLTexImage *otex = GetBaseLevel(i + _octave_min, DATA_ROT) + 2;
02478                 for(int j = 0; j < param._dog_level_num; j++, ftex++,  otex++, idx ++)
02479                 {
02480                         if(_levelFeatureNum[idx]<=0)continue;
02481                         sigma = param.GetLevelSigma(j+param._level_min+1);
02482                         //
02483                         ftex->AttachToFBO(0);
02484                         ftex->FitTexViewPort();
02485 
02486                         glActiveTexture(GL_TEXTURE0);
02487                         ftex->BindTex();
02488                         glActiveTexture(GL_TEXTURE1);
02489                         otex->BindTex();
02490 
02491                         ShaderMan::UseShaderSimpleOrientation(otex->GetTexID(), sigma, sigma_step);
02492                         ftex->DrawQuad();
02493                 }
02494         }
02495         GLTexImage::UnbindMultiTex(2);
02496 }
02497 
02498 void PyramidPacked::InitPyramid(int w, int h, int ds)
02499 {
02500         int wp, hp, toobig = 0;
02501         if(ds == 0)
02502         {
02503                 _down_sample_factor = 0;
02504                 if(GlobalUtil::_octave_min_default>=0)
02505                 {
02506                         wp = w >> GlobalUtil::_octave_min_default;
02507                         hp = h >> GlobalUtil::_octave_min_default;
02508                 }else 
02509                 {
02510                         wp = w << (-GlobalUtil::_octave_min_default);
02511                         hp = h << (-GlobalUtil::_octave_min_default);
02512                 }
02513                 _octave_min = _octave_min_default;
02514         }else
02515         {
02516                 //must use 0 as _octave_min; 
02517                 _octave_min = 0;
02518                 _down_sample_factor = ds;
02519                 w >>= ds;
02520                 h >>= ds;
02521                 wp = w;
02522                 hp = h; 
02523         }
02524 
02525         while(wp > GlobalUtil::_texMaxDim  || hp > GlobalUtil::_texMaxDim )
02526         {
02527                 _octave_min ++;
02528                 wp >>= 1;
02529                 hp >>= 1;
02530                 toobig = 1;
02531         }
02532         if(toobig && GlobalUtil::_verbose && _octave_min > 0)
02533         {
02534 
02535                 std::cout<< "**************************************************************\n"
02536                                         "Image larger than allowed dimension, data will be downsampled!\n"
02537                                         "use -maxd to change the settings\n"
02538                                         "***************************************************************\n";
02539         }
02540 
02541         if( wp == _pyramid_width && hp == _pyramid_height && _allocated )
02542         {
02543                 FitPyramid(wp, hp);
02544         }else if(GlobalUtil::_ForceTightPyramid || _allocated ==0)
02545         {
02546                 ResizePyramid(wp, hp);
02547         }
02548         else if( wp > _pyramid_width || hp > _pyramid_height )
02549         {
02550                 ResizePyramid(max(wp, _pyramid_width), max(hp, _pyramid_height));
02551                 if(wp < _pyramid_width || hp < _pyramid_height)  FitPyramid(wp, hp);
02552         }
02553         else
02554         {
02555                 //try use the pyramid allocated for large image on small input images
02556                 FitPyramid(wp, hp);
02557         }
02558 
02559         //select the initial smoothing filter according to the new _octave_min
02560         ShaderMan::SelectInitialSmoothingFilter(_octave_min + _down_sample_factor, param);
02561 }
02562 
02563 
02564 
02565 void PyramidPacked::FitPyramid(int w, int h)
02566 {
02567         //(w, h) <= (_pyramid_width, _pyramid_height);
02568 
02569         _pyramid_octave_first = 0;
02570         //
02571         _octave_num  = GlobalUtil::_octave_num_default;
02572 
02573         int _octave_num_max = GetRequiredOctaveNum(min(w, h));
02574 
02575         if(_octave_num < 1 || _octave_num > _octave_num_max) 
02576         {
02577                 _octave_num = _octave_num_max;
02578         }
02579 
02580 
02581         int pw = _pyramid_width>>1, ph = _pyramid_height>>1;
02582         while(_pyramid_octave_first + _octave_num < _pyramid_octave_num &&  
02583                 pw >= w && ph >= h)
02584         {
02585                 _pyramid_octave_first++;
02586                 pw >>= 1;
02587                 ph >>= 1;
02588         }
02589 
02590         for(int i = 0; i < _octave_num; i++)
02591         {
02592                 GLTexImage * tex = GetBaseLevel(i + _octave_min);
02593                 GLTexImage * dog = GetBaseLevel(i + _octave_min, DATA_DOG);
02594                 GLTexImage * grd = GetBaseLevel(i + _octave_min, DATA_GRAD);
02595                 GLTexImage * rot = GetBaseLevel(i + _octave_min, DATA_ROT);
02596                 GLTexImage * key = GetBaseLevel(i + _octave_min, DATA_KEYPOINT);
02597                 for(int j = param._level_min; j <= param._level_max; j++, tex++, dog++, grd++, rot++, key++)
02598                 {
02599                         tex->SetImageSize(w, h);
02600                         if(j == param._level_min) continue;
02601                         dog->SetImageSize(w, h);
02602                         grd->SetImageSize(w, h);
02603                         rot->SetImageSize(w, h);
02604                         if(j == param._level_min + 1 || j == param._level_max) continue;
02605                         key->SetImageSize(w, h);
02606                 }
02607                 w>>=1;
02608                 h>>=1;
02609         }
02610 }
02611 
02612 
02613 void PyramidPacked::ResizePyramid( int w,  int h)
02614 {
02615         //
02616         unsigned int totalkb = 0;
02617         int _octave_num_new, input_sz, i, j;
02618         //
02619 
02620         if(_pyramid_width == w && _pyramid_height == h && _allocated) return;
02621 
02622         if(w > GlobalUtil::_texMaxDim || h > GlobalUtil::_texMaxDim) return ;
02623 
02624         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<w<<"x"<<h<<endl;
02625         //first octave does not change
02626         _pyramid_octave_first = 0;
02627 
02628         
02629         //compute # of octaves
02630 
02631         input_sz = min(w,h) ;
02632         _pyramid_width =  w;
02633         _pyramid_height =  h;
02634 
02635         //reset to preset parameters
02636         _octave_num_new  = GlobalUtil::_octave_num_default;
02637 
02638         if(_octave_num_new < 1) _octave_num_new = GetRequiredOctaveNum(input_sz) ;
02639 
02640 
02641         if(_pyramid_octave_num != _octave_num_new)
02642         {
02643                 //destroy the original pyramid if the # of octave changes
02644                 if(_octave_num >0) 
02645                 {
02646                         DestroyPerLevelData();
02647                         DestroyPyramidData();
02648                 }
02649                 _pyramid_octave_num = _octave_num_new;
02650         }
02651 
02652         _octave_num = _pyramid_octave_num;
02653 
02654         int noct = _octave_num;
02655         int nlev = param._level_num;
02656 
02657         //      //initialize the pyramid
02658         if(_allPyramid==NULL)   _allPyramid = new GLTexPacked[ noct* nlev * DATA_NUM];
02659 
02660 
02661         GLTexPacked * gus = (GLTexPacked *) GetBaseLevel(_octave_min, DATA_GAUSSIAN);
02662         GLTexPacked * dog = (GLTexPacked *) GetBaseLevel(_octave_min, DATA_DOG);
02663         GLTexPacked * grd = (GLTexPacked *) GetBaseLevel(_octave_min, DATA_GRAD);
02664         GLTexPacked * rot = (GLTexPacked *) GetBaseLevel(_octave_min, DATA_ROT);
02665         GLTexPacked * key = (GLTexPacked *) GetBaseLevel(_octave_min, DATA_KEYPOINT);
02666 
02667 
02669 
02670         for(i = 0; i< noct; i++)
02671         {
02672                 for( j = 0; j< nlev; j++, gus++, dog++, grd++, rot++, key++)
02673                 {
02674                         gus->InitTexture(w, h);
02675                         if(j==0)continue;
02676                         dog->InitTexture(w, h);
02677                         grd->InitTexture(w, h, 0);
02678                         rot->InitTexture(w, h);
02679                         if(j<=1 || j >=nlev -1) continue;
02680                         key->InitTexture(w, h, 0);
02681                 }
02682                 int tsz = (gus -1)->GetTexPixelCount() * 16;
02683                 totalkb += ((nlev *5 -6)* tsz / 1024);
02684                 //several auxilary textures are not actually required
02685                 w>>=1;
02686                 h>>=1;
02687         }
02688 
02689         totalkb += ResizeFeatureStorage();
02690 
02691         _allocated = 1;
02692 
02693         if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<(totalkb/1024)<<"MB\n";
02694 
02695 }
02696 
02697 void PyramidPacked::DestroyPyramidData()
02698 {
02699         if(_allPyramid)
02700         {
02701                 delete [] _allPyramid;
02702                 _allPyramid = NULL;
02703         }
02704 }
02705 
02706 
02707 GLTexImage*  PyramidPacked::GetLevelTexture(int octave, int level, int dataName)
02708 {
02709         return _allPyramid+ (_pyramid_octave_first + octave - _octave_min) * param._level_num 
02710                 + param._level_num * _pyramid_octave_num * dataName
02711                 + (level - param._level_min);
02712 
02713 }
02714 
02715 GLTexImage*  PyramidPacked::GetLevelTexture(int octave, int level)
02716 {
02717         return _allPyramid+ (_pyramid_octave_first + octave - _octave_min) * param._level_num 
02718                 + (level - param._level_min);
02719 }
02720 
02721 //in the packed implementation( still in progress)
02722 // DATA_GAUSSIAN, DATA_DOG, DATA_GAD will be stored in different textures.
02723 
02724 GLTexImage*  PyramidPacked::GetBaseLevel(int octave, int dataName)
02725 {
02726         if(octave <_octave_min || octave > _octave_min + _octave_num) return NULL;
02727         int offset = (_pyramid_octave_first + octave - _octave_min) * param._level_num;
02728         int num = param._level_num * _pyramid_octave_num;
02729         return _allPyramid + num *dataName + offset;
02730 }
02731 
02732 
02733 void PyramidPacked::FitHistogramPyramid()
02734 {
02735         GLTexImage * tex, *htex;
02736         int hist_level_num = _hpLevelNum - _pyramid_octave_first; 
02737 
02738         tex = GetBaseLevel(_octave_min , DATA_KEYPOINT) + 2;
02739         htex = _histoPyramidTex + hist_level_num - 1;
02740         int w = (tex->GetImgWidth() + 2) >> 2;
02741         int h = (tex->GetImgHeight() + 2)>> 2;
02742 
02743 
02744         //4n+1 -> n; 4n+2,2, 3 -> n+1
02745         for(int k = 0; k <hist_level_num -1; k++, htex--)
02746         {
02747                 if(htex->GetImgHeight()!= h || htex->GetImgWidth() != w)
02748                 {       
02749                         htex->SetImageSize(w, h);
02750                         htex->ZeroHistoMargin();
02751                 }
02752 
02753                 w = (w + 1)>>1; h = (h + 1) >> 1;
02754         }
02755 }
02756 
 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