Statistics
| Branch: | Revision:

root / rgbdslam / external / siftgpu / src / SiftGPU / PyramidCL.cpp @ 9240aaa3

History | View | Annotate | Download (29.7 KB)

1
////////////////////////////////////////////////////////////////////////////
2
//        File:                PyramidCL.cpp
3
//        Author:                Changchang Wu
4
//        Description : implementation of the PyramidCL class.
5
//                                OpenCL-based implementation of SiftPyramid
6
//
7
//        Copyright (c) 2007 University of North Carolina at Chapel Hill
8
//        All Rights Reserved
9
//
10
//        Permission to use, copy, modify and distribute this software and its
11
//        documentation for educational, research and non-profit purposes, without
12
//        fee, and without a written agreement is hereby granted, provided that the
13
//        above copyright notice and the following paragraph appear in all copies.
14
//        
15
//        The University of North Carolina at Chapel Hill make no representations
16
//        about the suitability of this software for any purpose. It is provided
17
//        'as is' without express or implied warranty. 
18
//
19
//        Please send BUG REPORTS to ccwu@cs.unc.edu
20
//
21
////////////////////////////////////////////////////////////////////////////
22

    
23
#if defined(CL_SIFTGPU_ENABLED)
24

    
25

    
26
#include "GL/glew.h"
27
#include <CL/OpenCL.h>
28
#include <iostream>
29
#include <vector>
30
#include <algorithm>
31
#include <stdlib.h>
32
#include <math.h>
33
using namespace std;
34

    
35
#include "GlobalUtil.h"
36
#include "GLTexImage.h"
37
#include "CLTexImage.h" 
38
#include "SiftGPU.h"
39
#include "SiftPyramid.h"
40
#include "ProgramCL.h"
41
#include "PyramidCL.h"
42

    
43

    
44
#define USE_TIMING()                double t, t0, tt;
45
#define OCTAVE_START()                if(GlobalUtil::_timingO){        t = t0 = CLOCK();        cout<<"#"<<i+_down_sample_factor<<"\t";        }
46
#define LEVEL_FINISH()                if(GlobalUtil::_timingL){        _OpenCL->FinishCL();        tt = CLOCK();cout<<(tt-t)<<"\t";        t = CLOCK();}
47
#define OCTAVE_FINISH()                if(GlobalUtil::_timingO)cout<<"|\t"<<(CLOCK()-t0)<<endl;
48

    
49

    
50
PyramidCL::PyramidCL(SiftParam& sp) : SiftPyramid(sp)
51
{
52
        _allPyramid = NULL;
53
        _histoPyramidTex = NULL;
54
        _featureTex = NULL;
55
        _descriptorTex = NULL;
56
        _orientationTex = NULL;
57
    _bufferTEX = NULL;
58
    if(GlobalUtil::_usePackedTex)    _OpenCL = new ProgramBagCL();
59
    else                             _OpenCL = new ProgramBagCLN();
60
    _OpenCL->InitProgramBag(sp);
61
        _inputTex = new CLTexImage( _OpenCL->GetContextCL(),
62
                                _OpenCL->GetCommandQueue());
63
    /////////////////////////
64
    InitializeContext();
65
}
66

    
67
PyramidCL::~PyramidCL()
68
{
69
        DestroyPerLevelData();
70
        DestroySharedData();
71
        DestroyPyramidData();
72
    if(_OpenCL)   delete _OpenCL;
73
        if(_inputTex) delete _inputTex;
74
    if(_bufferTEX) delete _bufferTEX;
75
}
76

    
77
void PyramidCL::InitializeContext()
78
{
79
    GlobalUtil::InitGLParam(1);
80
}
81

    
82
void PyramidCL::InitPyramid(int w, int h, int ds)
83
{
84
        int wp, hp, toobig = 0;
85
        if(ds == 0)
86
        {
87
                _down_sample_factor = 0;
88
                if(GlobalUtil::_octave_min_default>=0)
89
                {
90
                        wp = w >> _octave_min_default;
91
                        hp = h >> _octave_min_default;
92
                }else
93
                {
94
                        //can't upsample by more than 8
95
                        _octave_min_default = max(-3, _octave_min_default);
96
                        //
97
                        wp = w << (-_octave_min_default);
98
                        hp = h << (-_octave_min_default);
99
                }
100
                _octave_min = _octave_min_default;
101
        }else
102
        {
103
                //must use 0 as _octave_min; 
104
                _octave_min = 0;
105
                _down_sample_factor = ds;
106
                w >>= ds;
107
                h >>= ds;
108
                wp = w;
109
                hp = h; 
110
        }
111

    
112
        while(wp > GlobalUtil::_texMaxDim  || hp > GlobalUtil::_texMaxDim )
113
        {
114
                _octave_min ++;
115
                wp >>= 1;
116
                hp >>= 1;
117
                toobig = 1;
118
        }
119
        if(toobig && GlobalUtil::_verbose && _octave_min > 0)
120
        {
121
                std::cout<< "**************************************************************\n"
122
                                        "Image larger than allowed dimension, data will be downsampled!\n"
123
                                        "use -maxd to change the settings\n"
124
                                        "***************************************************************\n";
125
        }
126

    
127
        if( wp == _pyramid_width && hp == _pyramid_height && _allocated )
128
        {
129
                FitPyramid(wp, hp);
130
        }else if(GlobalUtil::_ForceTightPyramid || _allocated ==0)
131
        {
132
                ResizePyramid(wp, hp);
133
        }
134
        else if( wp > _pyramid_width || hp > _pyramid_height )
135
        {
136
                ResizePyramid(max(wp, _pyramid_width), max(hp, _pyramid_height));
137
                if(wp < _pyramid_width || hp < _pyramid_height)  FitPyramid(wp, hp);
138
        }
139
        else
140
        {
141
                //try use the pyramid allocated for large image on small input images
142
                FitPyramid(wp, hp);
143
        }
144

    
145
    _OpenCL->SelectInitialSmoothingFilter(_octave_min + _down_sample_factor, param);
146
}
147

    
148
void PyramidCL::ResizePyramid(int w, int h)
149
{
150
        //
151
        unsigned int totalkb = 0;
152
        int _octave_num_new, input_sz, i, j;
153
        //
154

    
155
        if(_pyramid_width == w && _pyramid_height == h && _allocated) return;
156

    
157
        if(w > GlobalUtil::_texMaxDim || h > GlobalUtil::_texMaxDim) return ;
158

    
159
        if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<w<<"x"<<h<<endl;
160
        //first octave does not change
161
        _pyramid_octave_first = 0;
162

    
163
        
164
        //compute # of octaves
165
        input_sz = min(w,h) ;
166
        _pyramid_width =  w;
167
        _pyramid_height =  h;
168

    
169
        //reset to preset parameters
170
        _octave_num_new  = GlobalUtil::_octave_num_default;
171

    
172
        if(_octave_num_new < 1) _octave_num_new = GetRequiredOctaveNum(input_sz) ;
173

    
174
        if(_pyramid_octave_num != _octave_num_new)
175
        {
176
                //destroy the original pyramid if the # of octave changes
177
                if(_octave_num >0) 
178
                {
179
                        DestroyPerLevelData();
180
                        DestroyPyramidData();
181
                }
182
                _pyramid_octave_num = _octave_num_new;
183
        }
184

    
185
        _octave_num = _pyramid_octave_num;
186

    
187
        int noct = _octave_num;
188
        int nlev = param._level_num;
189
    int texNum = noct* nlev * DATA_NUM;
190

    
191
        //        //initialize the pyramid
192
        if(_allPyramid==NULL)
193
    {
194
        _allPyramid = new CLTexImage[ texNum];
195
        cl_context       context  = _OpenCL->GetContextCL();
196
        cl_command_queue queue    = _OpenCL->GetCommandQueue();
197
        for(i = 0; i < texNum; ++i)   _allPyramid[i].SetContext(context, queue);
198
    }
199

    
200

    
201

    
202
        CLTexImage * gus =  GetBaseLevel(_octave_min, DATA_GAUSSIAN);
203
        CLTexImage * dog =  GetBaseLevel(_octave_min, DATA_DOG);
204
        CLTexImage * grd =  GetBaseLevel(_octave_min, DATA_GRAD);
205
    CLTexImage * rot =  GetBaseLevel(_octave_min, DATA_ROT);
206
        CLTexImage * key =  GetBaseLevel(_octave_min, DATA_KEYPOINT);
207

    
208
        ////////////there could be "out of memory" happening during the allocation
209

    
210

    
211

    
212
        for(i = 0; i< noct; i++)
213
        {
214
                for( j = 0; j< nlev; j++, gus++, dog++, grd++, rot++, key++)
215
                {
216
            gus->InitPackedTex(w, h, GlobalUtil::_usePackedTex);
217
                        if(j==0)continue;
218
                        dog->InitPackedTex(w, h, GlobalUtil::_usePackedTex);
219
            if(j < 1 + param._dog_level_num)
220
                        {
221
                            grd->InitPackedTex(w, h, GlobalUtil::_usePackedTex);
222
                            rot->InitPackedTex(w, h, GlobalUtil::_usePackedTex);
223
            }
224
                        if(j > 1 && j < nlev -1) key->InitPackedTex(w, h, GlobalUtil::_usePackedTex);
225
                }
226
        ////////////////////////////////////////
227
                int tsz = (gus -1)->GetTexPixelCount() * 16;
228
                totalkb += ((nlev *5 -6)* tsz / 1024);
229
                //several auxilary textures are not actually required
230
                w>>=1;
231
                h>>=1;
232
        }
233

    
234
        totalkb += ResizeFeatureStorage();
235

    
236
        _allocated = 1;
237

    
238
        if(GlobalUtil::_verbose && GlobalUtil::_timingS) std::cout<<"[Allocate Pyramid]:\t" <<(totalkb/1024)<<"MB\n";
239

    
240
}
241

    
242
void PyramidCL::FitPyramid(int w, int h)
243
{
244
        _pyramid_octave_first = 0;
245
        //
246
        _octave_num  = GlobalUtil::_octave_num_default;
247

    
248
        int _octave_num_max = GetRequiredOctaveNum(min(w, h));
249

    
250
        if(_octave_num < 1 || _octave_num > _octave_num_max) 
251
        {
252
                _octave_num = _octave_num_max;
253
        }
254

    
255

    
256
        int pw = _pyramid_width>>1, ph = _pyramid_height>>1;
257
        while(_pyramid_octave_first + _octave_num < _pyramid_octave_num &&  
258
                pw >= w && ph >= h)
259
        {
260
                _pyramid_octave_first++;
261
                pw >>= 1;
262
                ph >>= 1;
263
        }
264

    
265
        //////////////////
266
        for(int i = 0; i < _octave_num; i++)
267
        {
268
                CLTexImage * tex = GetBaseLevel(i + _octave_min);
269
                CLTexImage * dog = GetBaseLevel(i + _octave_min, DATA_DOG);
270
                CLTexImage * grd = GetBaseLevel(i + _octave_min, DATA_GRAD);
271
                CLTexImage * rot = GetBaseLevel(i + _octave_min, DATA_ROT);
272
                CLTexImage * key = GetBaseLevel(i + _octave_min, DATA_KEYPOINT);
273
                for(int j = param._level_min; j <= param._level_max; j++, tex++, dog++, grd++, rot++, key++)
274
                {
275
                        tex->SetPackedSize(w, h, GlobalUtil::_usePackedTex);
276
                        if(j == param._level_min) continue;
277
                        dog->SetPackedSize(w, h, GlobalUtil::_usePackedTex);
278
            if(j < param._level_max - 1)
279
            {
280
                            grd->SetPackedSize(w, h, GlobalUtil::_usePackedTex);
281
                            rot->SetPackedSize(w, h, GlobalUtil::_usePackedTex);
282
            }
283
                        if(j > param._level_min + 1 &&  j < param._level_max) key->SetPackedSize(w, h, GlobalUtil::_usePackedTex);
284
                }
285
                w>>=1;
286
                h>>=1;
287
        }
288
}
289

    
290

    
291
void PyramidCL::SetLevelFeatureNum(int idx, int fcount)
292
{
293
    _featureTex[idx].InitBufferTex(fcount, 1, 4);
294
        _levelFeatureNum[idx] = fcount;
295
}
296

    
297
int PyramidCL::ResizeFeatureStorage()
298
{
299
        int totalkb = 0;
300
        if(_levelFeatureNum==NULL)        _levelFeatureNum = new int[_octave_num * param._dog_level_num];
301
        std::fill(_levelFeatureNum, _levelFeatureNum+_octave_num * param._dog_level_num, 0); 
302

    
303
    cl_context       context  = _OpenCL->GetContextCL();
304
    cl_command_queue queue    = _OpenCL->GetCommandQueue();
305
        int wmax = GetBaseLevel(_octave_min)->GetImgWidth() * 2;
306
        int hmax = GetBaseLevel(_octave_min)->GetImgHeight() * 2;
307
        int whmax = max(wmax, hmax);
308
        int w,  i;
309

    
310
        //
311
        int num = (int)ceil(log(double(whmax))/log(4.0));
312

    
313
        if( _hpLevelNum != num)
314
        {
315
                _hpLevelNum = num;
316
                if(_histoPyramidTex ) delete [] _histoPyramidTex;
317
                _histoPyramidTex = new CLTexImage[_hpLevelNum];
318
        for(i = 0; i < _hpLevelNum; ++i) _histoPyramidTex[i].SetContext(context, queue);
319
        }
320

    
321
        for(i = 0, w = 1; i < _hpLevelNum; i++)
322
        {
323
        _histoPyramidTex[i].InitBufferTex(w, whmax, 4);
324
                w<<=2;
325
        }
326

    
327
        // (4 ^ (_hpLevelNum) -1 / 3) pixels
328
        totalkb += (((1 << (2 * _hpLevelNum)) -1) / 3 * 16 / 1024);
329

    
330
        //initialize the feature texture
331
        int idx = 0, n = _octave_num * param._dog_level_num;
332
        if(_featureTex==NULL)        
333
    {
334
        _featureTex = new CLTexImage[n];
335
        for(i = 0; i <n; ++i) _featureTex[i].SetContext(context, queue);
336
    }
337
        if(GlobalUtil::_MaxOrientation >1 && GlobalUtil::_OrientationPack2==0 && _orientationTex== NULL)        
338
    {
339
        _orientationTex = new CLTexImage[n];
340
        for(i = 0; i < n; ++i) _orientationTex[i].SetContext(context, queue);
341
    }
342

    
343

    
344
        for(i = 0; i < _octave_num; i++)
345
        {
346
                CLTexImage * tex = GetBaseLevel(i+_octave_min);
347
                int fmax = int(4 * tex->GetTexWidth() * tex->GetTexHeight()*GlobalUtil::_MaxFeaturePercent);
348
                //
349
                if(fmax > GlobalUtil::_MaxLevelFeatureNum) fmax = GlobalUtil::_MaxLevelFeatureNum;
350
                else if(fmax < 32) fmax = 32;        //give it at least a space of 32 feature
351

    
352
                for(int j = 0; j < param._dog_level_num; j++, idx++)
353
                {
354
                        _featureTex[idx].InitBufferTex(fmax, 1, 4);
355
                        totalkb += fmax * 16 /1024;
356
                        //
357
                        if(GlobalUtil::_MaxOrientation>1 && GlobalUtil::_OrientationPack2 == 0)
358
                        {
359
                                _orientationTex[idx].InitBufferTex(fmax, 1, 4);
360
                                totalkb += fmax * 16 /1024;
361
                        }
362
                }
363
        }
364

    
365
        //this just need be initialized once
366
        if(_descriptorTex==NULL)
367
        {
368
                //initialize feature texture pyramid
369
                int fmax = _featureTex->GetImgWidth();
370
                _descriptorTex = new CLTexImage(context, queue);
371
                totalkb += ( fmax /2);
372
                _descriptorTex->InitBufferTex(fmax *128, 1, 1);
373
        }else
374
        {
375
                totalkb +=  _descriptorTex->GetDataSize()/1024;
376
        }
377
        return totalkb;
378
}
379

    
380
void PyramidCL::GetFeatureDescriptors() 
381
{
382
        //descriptors...
383
        /*float* pd =  &_descriptor_buffer[0];
384
        vector<float> descriptor_buffer2;
385

386
        //use another buffer if we need to re-order the descriptors
387
        if(_keypoint_index.size() > 0)
388
        {
389
                descriptor_buffer2.resize(_descriptor_buffer.size());
390
                pd = &descriptor_buffer2[0];
391
        }
392

393
        CLTexImage * got, * ftex= _featureTex;
394
        for(int i = 0, idx = 0; i < _octave_num; i++)
395
        {
396
                got = GetBaseLevel(i + _octave_min, DATA_GRAD) + 1;
397
                for(int j = 0; j < param._dog_level_num; j++, ftex++, idx++, got++)
398
                {
399
                        if(_levelFeatureNum[idx]==0) continue;
400
                        ProgramCL::ComputeDescriptor(ftex, got, _descriptorTex);//process
401
                        _descriptorTex->CopyToHost(pd); //readback descriptor
402
                        pd += 128*_levelFeatureNum[idx];
403
                }
404
        }
405

406
        if(GlobalUtil::_timingS) _OpenCL->FinishCL();
407

408
        if(_keypoint_index.size() > 0)
409
        {
410
            //put the descriptor back to the original order for keypoint list.
411
                for(int i = 0; i < _featureNum; ++i)
412
                {
413
                        int index = _keypoint_index[i];
414
                        memcpy(&_descriptor_buffer[index*128], &descriptor_buffer2[i*128], 128 * sizeof(float));
415
                }
416
        }*/ 
417
}
418

    
419
void PyramidCL::GenerateFeatureListTex() 
420
{
421

    
422
        vector<float> list;
423
        int idx = 0;
424
        const double twopi = 2.0*3.14159265358979323846;
425
        float sigma_half_step = powf(2.0f, 0.5f / param._dog_level_num);
426
        float octave_sigma = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
427
        float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f; 
428
        if(_down_sample_factor>0) octave_sigma *= float(1<<_down_sample_factor); 
429

    
430
        _keypoint_index.resize(0); // should already be 0
431
        for(int i = 0; i < _octave_num; i++, octave_sigma*= 2.0f)
432
        {
433
                for(int j = 0; j < param._dog_level_num; j++, idx++)
434
                {
435
                        list.resize(0);
436
                        float level_sigma = param.GetLevelSigma(j + param._level_min + 1) * octave_sigma;
437
                        float sigma_min = level_sigma / sigma_half_step;
438
                        float sigma_max = level_sigma * sigma_half_step;
439
                        int fcount = 0 ;
440
                        for(int k = 0; k < _featureNum; k++)
441
                        {
442
                                float * key = &_keypoint_buffer[k*4];
443
                                if(   (key[2] >= sigma_min && key[2] < sigma_max)
444
                                        ||(key[2] < sigma_min && i ==0 && j == 0)
445
                                        ||(key[2] > sigma_max && i == _octave_num -1 && j == param._dog_level_num - 1))
446
                                {
447
                                        //add this keypoint to the list
448
                                        list.push_back((key[0] - offset) / octave_sigma + 0.5f);
449
                                        list.push_back((key[1] - offset) / octave_sigma + 0.5f);
450
                                        list.push_back(key[2] / octave_sigma);
451
                                        list.push_back((float)fmod(twopi-key[3], twopi));
452
                                        fcount ++;
453
                                        //save the index of keypoints
454
                                        _keypoint_index.push_back(k);
455
                                }
456

    
457
                        }
458

    
459
                        _levelFeatureNum[idx] = fcount;
460
                        if(fcount==0)continue;
461
                        CLTexImage * ftex = _featureTex+idx;
462

    
463
                        SetLevelFeatureNum(idx, fcount);
464
                        ftex->CopyFromHost(&list[0]);
465
                }
466
        }
467

    
468
        if(GlobalUtil::_verbose)
469
        {
470
                std::cout<<"#Features:\t"<<_featureNum<<"\n";
471
        }
472

    
473
}
474

    
475
void PyramidCL::ReshapeFeatureListCPU() 
476
{
477
        int i, szmax =0, sz;
478
        int n = param._dog_level_num*_octave_num;
479
        for( i = 0; i < n; i++) 
480
        {
481
                sz = _levelFeatureNum[i];
482
                if(sz > szmax ) szmax = sz;
483
        }
484
        float * buffer = new float[szmax*16];
485
        float * buffer1 = buffer;
486
        float * buffer2 = buffer + szmax*4;
487

    
488

    
489

    
490
        _featureNum = 0;
491

    
492
#ifdef NO_DUPLICATE_DOWNLOAD
493
        const double twopi = 2.0*3.14159265358979323846;
494
        _keypoint_buffer.resize(0);
495
        float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
496
        if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
497
        float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
498
#endif
499

    
500

    
501
        for(i = 0; i < n; i++)
502
        {
503
                if(_levelFeatureNum[i]==0)continue;
504

    
505
                _featureTex[i].CopyToHost(buffer1);
506
                
507
                int fcount =0;
508
                float * src = buffer1;
509
                float * des = buffer2;
510
                const static double factor  = 2.0*3.14159265358979323846/65535.0;
511
                for(int j = 0; j < _levelFeatureNum[i]; j++, src+=4)
512
                {
513
                        unsigned short * orientations = (unsigned short*) (&src[3]);
514
                        if(orientations[0] != 65535)
515
                        {
516
                                des[0] = src[0];
517
                                des[1] = src[1];
518
                                des[2] = src[2];
519
                                des[3] = float( factor* orientations[0]);
520
                                fcount++;
521
                                des += 4;
522
                                if(orientations[1] != 65535 && orientations[1] != orientations[0])
523
                                {
524
                                        des[0] = src[0];
525
                                        des[1] = src[1];
526
                                        des[2] = src[2];
527
                                        des[3] = float(factor* orientations[1]);        
528
                                        fcount++;
529
                                        des += 4;
530
                                }
531
                        }
532
                }
533
                //texture size
534
                SetLevelFeatureNum(i, fcount);
535
                _featureTex[i].CopyFromHost(buffer2);
536
                
537
                if(fcount == 0) continue;
538

    
539
#ifdef NO_DUPLICATE_DOWNLOAD
540
                float oss = os * (1 << (i / param._dog_level_num));
541
                _keypoint_buffer.resize((_featureNum + fcount) * 4);
542
                float* ds = &_keypoint_buffer[_featureNum * 4];
543
                float* fs = buffer2;
544
                for(int k = 0;  k < fcount; k++, ds+=4, fs+=4)
545
                {
546
                        ds[0] = oss*(fs[0]-0.5f) + offset;        //x
547
                        ds[1] = oss*(fs[1]-0.5f) + offset;        //y
548
                        ds[2] = oss*fs[2];  //scale
549
                        ds[3] = (float)fmod(twopi-fs[3], twopi);        //orientation, mirrored
550
                }
551
#endif
552
                _featureNum += fcount;
553
        }
554
        delete[] buffer;
555
        if(GlobalUtil::_verbose)
556
        {
557
                std::cout<<"#Features MO:\t"<<_featureNum<<endl;
558
        }
559
}
560

    
561
void PyramidCL::GenerateFeatureDisplayVBO() 
562
{
563
        //it is weried that this part is very slow.
564
        //use a big VBO to save all the SIFT box vertices
565
        /*int nvbo = _octave_num * param._dog_level_num;
566
        if(_featureDisplayVBO==NULL)
567
        {
568
                //initialize the vbos
569
                _featureDisplayVBO = new GLuint[nvbo];
570
                _featurePointVBO = new GLuint[nvbo];
571
                glGenBuffers(nvbo, _featureDisplayVBO);        
572
                glGenBuffers(nvbo, _featurePointVBO);
573
        }
574
        for(int i = 0; i < nvbo; i++)
575
        {
576
                if(_levelFeatureNum[i]<=0)continue;
577
                CLTexImage * ftex  = _featureTex + i;
578
                CLTexImage texPBO1( _levelFeatureNum[i]* 10, 1, 4, _featureDisplayVBO[i]);
579
                CLTexImage texPBO2(_levelFeatureNum[i], 1, 4, _featurePointVBO[i]);
580
                _OpenCL->DisplayKeyBox(ftex, &texPBO1);
581
                _OpenCL->DisplayKeyPoint(ftex, &texPBO2);        
582
        }*/
583
}
584

    
585
void PyramidCL::DestroySharedData() 
586
{
587
        //histogram reduction
588
        if(_histoPyramidTex)
589
        {
590
                delete[]        _histoPyramidTex;
591
                _hpLevelNum = 0;
592
                _histoPyramidTex = NULL;
593
        }
594
        //descriptor storage shared by all levels
595
        if(_descriptorTex)
596
        {
597
                delete _descriptorTex;
598
                _descriptorTex = NULL;
599
        }
600
        //cpu reduction buffer.
601
        if(_histo_buffer)
602
        {
603
                delete[] _histo_buffer;
604
                _histo_buffer = 0;
605
        }
606
}
607

    
608
void PyramidCL::DestroyPerLevelData() 
609
{
610
        //integers vector to store the feature numbers.
611
        if(_levelFeatureNum)
612
        {
613
                delete [] _levelFeatureNum;
614
                _levelFeatureNum = NULL;
615
        }
616
        //texture used to store features
617
        if(        _featureTex)
618
        {
619
                delete [] _featureTex;
620
                _featureTex =        NULL;
621
        }
622
        //texture used for multi-orientation 
623
        if(_orientationTex)
624
        {
625
                delete [] _orientationTex;
626
                _orientationTex = NULL;
627
        }
628
        int no = _octave_num* param._dog_level_num;
629

    
630
        //two sets of vbos used to display the features
631
        if(_featureDisplayVBO)
632
        {
633
                glDeleteBuffers(no, _featureDisplayVBO);
634
                delete [] _featureDisplayVBO;
635
                _featureDisplayVBO = NULL;
636
        }
637
        if( _featurePointVBO)
638
        {
639
                glDeleteBuffers(no, _featurePointVBO);
640
                delete [] _featurePointVBO;
641
                _featurePointVBO = NULL;
642
        }
643
}
644

    
645
void PyramidCL::DestroyPyramidData()
646
{
647
        if(_allPyramid)
648
        {
649
                delete [] _allPyramid;
650
                _allPyramid = NULL;
651
        }
652
}
653

    
654
void PyramidCL::DownloadKeypoints() 
655
{
656
        const double twopi = 2.0*3.14159265358979323846;
657
        int idx = 0;
658
        float * buffer = &_keypoint_buffer[0];
659
        vector<float> keypoint_buffer2;
660
        //use a different keypoint buffer when processing with an exisint features list
661
        //without orientation information. 
662
        if(_keypoint_index.size() > 0)
663
        {
664
                keypoint_buffer2.resize(_keypoint_buffer.size());
665
                buffer = &keypoint_buffer2[0];
666
        }
667
        float * p = buffer, *ps;
668
        CLTexImage * ftex = _featureTex;
669
        /////////////////////
670
        float os = _octave_min>=0? float(1<<_octave_min): 1.0f/(1<<(-_octave_min));
671
        if(_down_sample_factor>0) os *= float(1<<_down_sample_factor); 
672
        float offset = GlobalUtil::_LoweOrigin? 0 : 0.5f;
673
        /////////////////////
674
        for(int i = 0; i < _octave_num; i++, os *= 2.0f)
675
        {
676
        
677
                for(int j = 0; j  < param._dog_level_num; j++, idx++, ftex++)
678
                {
679

    
680
                        if(_levelFeatureNum[idx]>0)
681
                        {        
682
                                ftex->CopyToHost(ps = p);
683
                                for(int k = 0;  k < _levelFeatureNum[idx]; k++, ps+=4)
684
                                {
685
                                        ps[0] = os*(ps[0]-0.5f) + offset;        //x
686
                                        ps[1] = os*(ps[1]-0.5f) + offset;        //y
687
                                        ps[2] = os*ps[2]; 
688
                                        ps[3] = (float)fmod(twopi-ps[3], twopi);        //orientation, mirrored
689
                                }
690
                                p+= 4* _levelFeatureNum[idx];
691
                        }
692
                }
693
        }
694

    
695
        //put the feature into their original order for existing keypoint 
696
        if(_keypoint_index.size() > 0)
697
        {
698
                for(int i = 0; i < _featureNum; ++i)
699
                {
700
                        int index = _keypoint_index[i];
701
                        memcpy(&_keypoint_buffer[index*4], &keypoint_buffer2[i*4], 4 * sizeof(float));
702
                }
703
        }
704
}
705

    
706
void PyramidCL::GenerateFeatureListCPU()
707
{
708
        //no cpu version provided
709
        GenerateFeatureList();
710
}
711

    
712
void PyramidCL::GenerateFeatureList(int i, int j, int reduction_count, vector<int>& hbuffer)
713
{
714
    /*int fcount = 0, idx = i * param._dog_level_num  + j;
715
        int hist_level_num = _hpLevelNum - _pyramid_octave_first /2; 
716
        int ii, k, len; 
717

718
        CLTexImage * htex, * ftex, * tex, *got;
719
        ftex = _featureTex + idx;
720
        htex = _histoPyramidTex + hist_level_num -1;
721
        tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2 + j;
722
        got = GetBaseLevel(_octave_min + i, DATA_GRAD) + 2 + j;
723

724
        _OpenCL->InitHistogram(tex, htex);
725

726
        for(k = 0; k < reduction_count - 1; k++, htex--)
727
        {
728
                ProgramCL::ReduceHistogram(htex, htex -1);        
729
        }
730
        
731
        //htex has the row reduction result
732
        len = htex->GetImgHeight() * 4;
733
        hbuffer.resize(len);
734
        _OpenCL->FinishCL();
735
        htex->CopyToHost(&hbuffer[0]);
736
        //
737
        for(ii = 0; ii < len; ++ii)                fcount += hbuffer[ii];
738
        SetLevelFeatureNum(idx, fcount);
739
        
740
    //build the feature list
741
        if(fcount > 0)
742
        {
743
                _featureNum += fcount;
744
                _keypoint_buffer.resize(fcount * 4);
745
                //vector<int> ikbuf(fcount*4);
746
                int* ibuf = (int*) (&_keypoint_buffer[0]);
747

748
                for(ii = 0; ii < len; ++ii)
749
                {
750
                        int x = ii%4, y = ii / 4;
751
                        for(int jj = 0 ; jj < hbuffer[ii]; ++jj, ibuf+=4)
752
                        {
753
                                ibuf[0] = x; ibuf[1] = y; ibuf[2] = jj; ibuf[3] = 0;
754
                        }
755
                }
756
                _featureTex[idx].CopyFromHost(&_keypoint_buffer[0]);
757
        
758
                ////////////////////////////////////////////
759
                ProgramCL::GenerateList(_featureTex + idx, ++htex);
760
                for(k = 2; k < reduction_count; k++)
761
                {
762
                        ProgramCL::GenerateList(_featureTex + idx, ++htex);
763
                }
764
        }*/
765
}
766

    
767
void PyramidCL::GenerateFeatureList()
768
{
769
        /*double t1, t2; 
770
        int ocount = 0, reduction_count;
771
    int reverse = (GlobalUtil::_TruncateMethod == 1);
772

773
        vector<int> hbuffer;
774
        _featureNum = 0;
775

776
        //for(int i = 0, idx = 0; i < _octave_num; i++)
777
    FOR_EACH_OCTAVE(i, reverse)
778
        {
779
        CLTexImage* tex = GetBaseLevel(_octave_min + i, DATA_KEYPOINT) + 2;
780
                reduction_count = FitHistogramPyramid(tex);
781

782
                if(GlobalUtil::_timingO)
783
                {
784
                        t1 = CLOCK(); 
785
                        ocount = 0;
786
                        std::cout<<"#"<<i+_octave_min + _down_sample_factor<<":\t";
787
                }
788
                //for(int j = 0; j < param._dog_level_num; j++, idx++)
789
        FOR_EACH_LEVEL(j, reverse)
790
                {
791
            if(GlobalUtil::_TruncateMethod && GlobalUtil::_FeatureCountThreshold > 0 && _featureNum > GlobalUtil::_FeatureCountThreshold) continue;
792

793
                GenerateFeatureList(i, j, reduction_count, hbuffer);
794

795
                        /////////////////////////////
796
                        if(GlobalUtil::_timingO)
797
                        {
798
                int idx = i * param._dog_level_num + j;
799
                                ocount += _levelFeatureNum[idx];
800
                                std::cout<< _levelFeatureNum[idx] <<"\t";
801
                        }
802
                }
803
                if(GlobalUtil::_timingO)
804
                {        
805
                        t2 = CLOCK(); 
806
                        std::cout << "| \t" << int(ocount) << " :\t(" << (t2 - t1) << ")\n";
807
                }
808
        }
809
        /////
810
        CopyGradientTex();
811
        /////
812
        if(GlobalUtil::_timingS)_OpenCL->FinishCL();
813

814
        if(GlobalUtil::_verbose)
815
        {
816
                std::cout<<"#Features:\t"<<_featureNum<<"\n";
817
        }*/
818
}
819

    
820
GLTexImage* PyramidCL::GetLevelTexture(int octave, int level)
821
{
822
        return GetLevelTexture(octave, level, DATA_GAUSSIAN);
823
}
824

    
825
GLTexImage* PyramidCL::ConvertTexCL2GL(CLTexImage* tex, int dataName)
826
{
827
   
828
    if(_bufferTEX == NULL) _bufferTEX = new GLTexImage;
829

    
830
    ///////////////////////////////////////////
831
    int ratio = GlobalUtil::_usePackedTex ? 2 : 1; 
832
    int width  = tex->GetImgWidth() * ratio;
833
    int height = tex->GetImgHeight() * ratio; 
834
    int tw = max(width,  _bufferTEX->GetTexWidth());
835
    int th = max(height, _bufferTEX->GetTexHeight());
836
    _bufferTEX->InitTexture(tw, th, 1, GL_RGBA);
837
    _bufferTEX->SetImageSize(width, height); 
838

    
839
    //////////////////////////////////
840
    CLTexImage texCL(_OpenCL->GetContextCL(), _OpenCL->GetCommandQueue()); 
841
    texCL.InitTextureGL(*_bufferTEX, width, height, 4); 
842

    
843
        switch(dataName)
844
        {
845
        case DATA_GAUSSIAN: _OpenCL->UnpackImage(tex, &texCL); break;
846
        case DATA_DOG:_OpenCL->UnpackImageDOG(tex, &texCL); break;
847
        case DATA_GRAD:_OpenCL->UnpackImageGRD(tex, &texCL); break;
848
        case DATA_KEYPOINT:_OpenCL->UnpackImageKEY(tex,
849
        tex - param._level_num * _pyramid_octave_num, &texCL);break;
850
        default:
851
                        break;
852
        }
853

    
854

    
855
        return _bufferTEX;
856
}
857

    
858
GLTexImage* PyramidCL::GetLevelTexture(int octave, int level, int dataName) 
859
{
860
        CLTexImage* tex = GetBaseLevel(octave, dataName) + (level - param._level_min);
861
        return ConvertTexCL2GL(tex, dataName);
862
}
863

    
864
void PyramidCL::ConvertInputToCL(GLTexInput* input, CLTexImage* output)
865
{
866
        int ws = input->GetImgWidth(), hs = input->GetImgHeight();
867
        //copy the input image to pixel buffer object
868
    if(input->_pixel_data)
869
    {
870
        output->InitTexture(ws, hs, 1);
871
        output->CopyFromHost(input->_pixel_data); 
872
    }else /*if(input->_rgb_converted && input->CopyToPBO(_bufferPBO, ws, hs, GL_LUMINANCE))
873
    {
874
                output->InitTexture(ws, hs, 1);
875
        output->CopyFromPBO(ws, hs, _bufferPBO); 
876
    }else if(input->CopyToPBO(_bufferPBO, ws, hs))
877
        {
878
                CLTexImage texPBO(ws, hs, 4, _bufferPBO);
879
                output->InitTexture(ws, hs, 1);
880
                ProgramCL::ReduceToSingleChannel(output, &texPBO, !input->_rgb_converted);
881
        }else*/
882
        {
883
                std::cerr<< "Unable To Convert Intput\n";
884
        }
885
}
886

    
887
void PyramidCL::BuildPyramid(GLTexInput * input)
888
{
889

    
890
        USE_TIMING();
891

    
892
        int i, j;
893
        
894
        for ( i = _octave_min; i < _octave_min + _octave_num; i++)
895
        {
896

    
897
                CLTexImage *tex = GetBaseLevel(i);
898
                CLTexImage *buf = GetBaseLevel(i, DATA_DOG) +2;
899
        FilterCL ** filter  = _OpenCL->f_gaussian_step; 
900
                j = param._level_min + 1;
901

    
902
                OCTAVE_START();
903

    
904
                if( i == _octave_min )
905
                {        
906
            if(GlobalUtil::_usePackedTex)
907
            {
908
                            ConvertInputToCL(input, _inputTex);
909
                            if(i < 0)        _OpenCL->SampleImageU(tex, _inputTex, -i- 1);                        
910
                            else                _OpenCL->SampleImageD(tex, _inputTex, i + 1);
911
            }else
912
            {
913
                if(i == 0) ConvertInputToCL(input, tex);
914
                else
915
                {
916
                    ConvertInputToCL(input, _inputTex);
917
                    if(i < 0) _OpenCL->SampleImageU(tex, _inputTex, -i);
918
                    else      _OpenCL->SampleImageD(tex, _inputTex,  i);
919
                }
920
            }
921
            _OpenCL->FilterInitialImage(tex, buf);   
922
                }else
923
                {
924
                        _OpenCL->SampleImageD(tex, GetBaseLevel(i - 1) + param._level_ds - param._level_min); 
925
            _OpenCL->FilterSampledImage(tex, buf);
926
                }
927
                LEVEL_FINISH();
928
                for( ; j <=  param._level_max ; j++, tex++, filter++)
929
                {
930
                        // filtering
931
                        _OpenCL->FilterImage(*filter, tex + 1, tex, buf);
932
                        LEVEL_FINISH();
933
                }
934
                OCTAVE_FINISH();
935
        }
936
        if(GlobalUtil::_timingS) _OpenCL->FinishCL();
937
}
938

    
939
void PyramidCL::DetectKeypointsEX()
940
{
941
        int i, j;
942
        double t0, t, ts, t1, t2;
943
    
944
        if(GlobalUtil::_timingS && GlobalUtil::_verbose) ts = CLOCK();
945

    
946
        for(i = _octave_min; i < _octave_min + _octave_num; i++)
947
        {
948
                CLTexImage * gus = GetBaseLevel(i) + 1;
949
                CLTexImage * dog = GetBaseLevel(i, DATA_DOG) + 1;
950
                CLTexImage * grd = GetBaseLevel(i, DATA_GRAD) + 1;
951
        CLTexImage * rot = GetBaseLevel(i, DATA_ROT) + 1;
952
                //compute the gradient
953
                for(j = param._level_min +1; j <=  param._level_max ; j++, gus++, dog++, grd++, rot++)
954
                {
955
                        //input: gus and gus -1
956
                        //output: gradient, dog, orientation
957
                        _OpenCL->ComputeDOG(gus, gus - 1, dog, grd, rot);
958
                }
959
        }
960
        if(GlobalUtil::_timingS && GlobalUtil::_verbose)
961
        {
962
                _OpenCL->FinishCL();
963
                t1 = CLOCK();
964
        }
965
    //if(GlobalUtil::_timingS) _OpenCL->FinishCL();
966
    //if(!GlobalUtil::_usePackedTex) return; //not finished
967
    //return; 
968

    
969
        for ( i = _octave_min; i < _octave_min + _octave_num; i++)
970
        {
971
                if(GlobalUtil::_timingO)
972
                {
973
                        t0 = CLOCK();
974
                        std::cout<<"#"<<(i + _down_sample_factor)<<"\t";
975
                }
976
                CLTexImage * dog = GetBaseLevel(i, DATA_DOG) + 2;
977
                CLTexImage * key = GetBaseLevel(i, DATA_KEYPOINT) +2;
978

    
979

    
980
                for( j = param._level_min +2; j <  param._level_max ; j++, dog++, key++)
981
                {
982
                        if(GlobalUtil::_timingL)t = CLOCK();
983
                        //input, dog, dog + 1, dog -1
984
                        //output, key
985
                        _OpenCL->ComputeKEY(dog, key, param._dog_threshold, param._edge_threshold);
986
                        if(GlobalUtil::_timingL)
987
                        {
988
                                std::cout<<(CLOCK()-t)<<"\t";
989
                        }
990
                }
991
                if(GlobalUtil::_timingO)
992
                {
993
                        std::cout<<"|\t"<<(CLOCK()-t0)<<"\n";
994
                }
995
        }
996

    
997
        if(GlobalUtil::_timingS)
998
        {
999
                _OpenCL->FinishCL();
1000
                if(GlobalUtil::_verbose) 
1001
                {        
1002
                        t2 = CLOCK();
1003
                        std::cout        <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n"
1004
                                                <<"<Get Keypoints  >\t"<<(t2-t1)<<"\n";
1005
                }                                
1006
        }
1007
}
1008

    
1009
void PyramidCL::CopyGradientTex()
1010
{
1011
        /*double ts, t1;
1012

1013
        if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
1014

1015
        for(int i = 0, idx = 0; i < _octave_num; i++)
1016
        {
1017
                CLTexImage * got = GetBaseLevel(i + _octave_min, DATA_GRAD) +  1;
1018
                //compute the gradient
1019
                for(int j = 0; j <  param._dog_level_num ; j++, got++, idx++)
1020
                {
1021
                        if(_levelFeatureNum[idx] > 0)        got->CopyToTexture2D();
1022
                }
1023
        }
1024
        if(GlobalUtil::_timingS)
1025
        {
1026
                ProgramCL::FinishCLDA();
1027
                if(GlobalUtil::_verbose)
1028
                {
1029
                        t1 = CLOCK();
1030
                        std::cout        <<"<Copy Grad/Orientation>\t"<<(t1-ts)<<"\n";
1031
                }
1032
        }*/
1033
}
1034

    
1035
void PyramidCL::ComputeGradient() 
1036
{
1037

    
1038
        /*int i, j;
1039
        double ts, t1;
1040

1041
        if(GlobalUtil::_timingS && GlobalUtil::_verbose)ts = CLOCK();
1042

1043
        for(i = _octave_min; i < _octave_min + _octave_num; i++)
1044
        {
1045
                CLTexImage * gus = GetBaseLevel(i) +  1;
1046
                CLTexImage * dog = GetBaseLevel(i, DATA_DOG) +  1;
1047
                CLTexImage * got = GetBaseLevel(i, DATA_GRAD) +  1;
1048

1049
                //compute the gradient
1050
                for(j = 0; j <  param._dog_level_num ; j++, gus++, dog++, got++)
1051
                {
1052
                        ProgramCL::ComputeDOG(gus, dog, got);
1053
                }
1054
        }
1055
        if(GlobalUtil::_timingS)
1056
        {
1057
                ProgramCL::FinishCLDA();
1058
                if(GlobalUtil::_verbose)
1059
                {
1060
                        t1 = CLOCK();
1061
                        std::cout        <<"<Gradient, DOG  >\t"<<(t1-ts)<<"\n";
1062
                }
1063
        }*/
1064
}
1065

    
1066
int PyramidCL::FitHistogramPyramid(CLTexImage* tex)
1067
{
1068
        CLTexImage *htex;
1069
        int hist_level_num = _hpLevelNum - _pyramid_octave_first / 2; 
1070
        htex = _histoPyramidTex + hist_level_num - 1;
1071
        int w = (tex->GetImgWidth() + 2) >> 2;
1072
        int h = tex->GetImgHeight();
1073
        int count = 0; 
1074
        for(int k = 0; k < hist_level_num; k++, htex--)
1075
        {
1076
                //htex->SetImageSize(w, h);        
1077
                htex->InitTexture(w, h, 4); 
1078
                ++count;
1079
                if(w == 1)
1080
            break;
1081
                w = (w + 3)>>2; 
1082
        }
1083
        return count;
1084
}
1085

    
1086
void PyramidCL::GetFeatureOrientations() 
1087
{
1088

    
1089
/*
1090
        CLTexImage * ftex = _featureTex;
1091
        int * count         = _levelFeatureNum;
1092
        float sigma, sigma_step = powf(2.0f, 1.0f/param._dog_level_num);
1093

1094
        for(int i = 0; i < _octave_num; i++)
1095
        {
1096
                CLTexImage* got = GetBaseLevel(i + _octave_min, DATA_GRAD) + 1;
1097
                CLTexImage* key = GetBaseLevel(i + _octave_min, DATA_KEYPOINT) + 2;
1098

1099
                for(int j = 0; j < param._dog_level_num; j++, ftex++, count++, got++, key++)
1100
                {
1101
                        if(*count<=0)continue;
1102

1103
                        //if(ftex->GetImgWidth() < *count) ftex->InitTexture(*count, 1, 4);
1104

1105
                        sigma = param.GetLevelSigma(j+param._level_min+1);
1106

1107
                        ProgramCL::ComputeOrientation(ftex, got, key, sigma, sigma_step, _existing_keypoints);                
1108
                }
1109
        }
1110

1111
        if(GlobalUtil::_timingS)ProgramCL::FinishCL();
1112
        */
1113

    
1114

    
1115
}
1116

    
1117
void PyramidCL::GetSimplifiedOrientation() 
1118
{
1119
        //no simplified orientation
1120
        GetFeatureOrientations();
1121
}
1122

    
1123
CLTexImage* PyramidCL::GetBaseLevel(int octave, int dataName)
1124
{
1125
        if(octave <_octave_min || octave > _octave_min + _octave_num) return NULL;
1126
        int offset = (_pyramid_octave_first + octave - _octave_min) * param._level_num;
1127
        int num = param._level_num * _pyramid_octave_num;
1128
        return _allPyramid + num * dataName + offset;
1129
}
1130

    
1131
#endif
1132