pyramidgl.cpp

SiftGPU is an implementation of SIFT [1] for GPU. SiftGPU processes pixels parallely to build Gaussi

////////////////////////////////////////////////////////////////////////////
//	File:		PyramidGL.cpp
//	Author:		Changchang Wu
//	Description : implementation of PyramidGL/PyramidNaive/PyramidPackdc .
//
//
//	Copyright (c) 2007 University of North Carolina at Chapel Hill
//	All Rights Reserved
//
//	Permission to use, copy, modify and distribute this software and its
//	documentation for educational, research and non-profit purposes, without
//	fee, and without a written agreement is hereby granted, provided that the
//	above copyright notice and the following paragraph appear in all copies.
//	
//	The University of North Carolina at Chapel Hill make no representations
//	about the suitability of this software for any purpose. It is provided
//	'as is' without express or implied warranty. 
//
//	Please send BUG REPORTS to [e-mail]
//
////////////////////////////////////////////////////////////////////////////

#include "GL/glew.h"
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <fstream>
#include <math.h>
using namespace std;

#include "GlobalUtil.h"
#include "GLTexImage.h"
#include "SiftGPU.h"
#include "ShaderMan.h"
#include "SiftPyramid.h"
#include "PyramidGL.h"
#include "FrameBufferObject.h"


#if defined(__SSE__) || _MSC_VER > 1200
#define USE_SSE_FOR_SIFTGPU
#include <xmmintrin.h>
#else
#pragma message( "SSE optimization off!\n" )
#endif


#define USE_TIMING()		double t, t0, tt;
#define OCTAVE_START()		if(GlobalUtil::_timingO){	t = t0 = CLOCK();	cout<<"#"<<i+_down_sample_factor<<"\t";	}
#define LEVEL_FINISH()		if(GlobalUtil::_timingL){	glFinish();	tt = CLOCK();cout<<(tt-t)<<"\t";	t = CLOCK();}
#define OCTAVE_FINISH()		if(GlobalUtil::_timingO)cout<<"|\t"<<(CLOCK()-t0)<<endl;


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
PyramidNaive::PyramidNaive(SiftParam& sp): PyramidGL(sp)
{
	_texPyramid = NULL;
	_auxPyramid = NULL;
}

PyramidNaive::~PyramidNaive()
{
	DestroyPyramidData();
}

//align must be 2^i
void PyramidGL::	GetAlignedStorageSize(int num, int align,  int &fw, int &fh)
{
	if(num <=0)
	{
		fw = fh = 0;
	}else if(num < align*align)
	{
		fw = align;
		fh = (int)ceil(double(num) / fw);
	}else if(GlobalUtil::_NarrowFeatureTex)
	{
		double dn = double(num);
		int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/align);	
		fw = align * nb;	
		fh = (int)ceil(dn /fw);/**/
	}else
	{
		double dn = double(num);
		int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/align);
		fh = align * nb;
		if(nb <=1)
		{
			fw = (int)ceil(dn / fh);
			//align this dimension to blocksize
			fw = ((int) ceil(double(fw) /align))*align;
		}else
		{
			fw = GlobalUtil::_texMaxDim;
		}

	}


}

void PyramidGL::GetTextureStorageSize(int num, int &fw, int& fh)
{
	if(num <=0)
	{
		fw = fh = 0;
	}else if(num <= GlobalUtil::_FeatureTexBlock)
	{
		fw = num;
		fh = 1;
	}else if(GlobalUtil::_NarrowFeatureTex)
	{
		double dn = double(num);
		int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/GlobalUtil::_FeatureTexBlock);	
		fw = GlobalUtil::_FeatureTexBlock * nb;	
		fh = (int)ceil(dn /fw);/**/
	}else
	{
		double dn = double(num);
		int nb = (int) ceil(dn/GlobalUtil::_texMaxDim/GlobalUtil::_FeatureTexBlock);
		fh = GlobalUtil::_FeatureTexBlock * nb;
		if(nb <=1)
		{
			fw = (int)ceil(dn / fh);

			//align this dimension to blocksize

			//
			if( fw < fh)
			{
				int temp = fh;
				fh = fw;
				fw = temp;
			}
		}else
		{
			fw = GlobalUtil::_texMaxDim;
		}
	}
}

void PyramidNaive::DestroyPyramidData()
{
	if(_texPyramid)
	{
		delete [] _texPyramid;
		_texPyramid = NULL;
	}
	if(_auxPyramid)
	{
		delete [] _auxPyramid;  
		_auxPyramid = NULL;
	}
}

PyramidGL::~PyramidGL()
{
	DestroyPerLevelData();
	DestroySharedData();
}

void PyramidGL::DestroyPerLevelData()
{
	//integers vector to store the feature numbers.
	if(_levelFeatureNum)
	{
		delete [] _levelFeatureNum;
		_levelFeatureNum = NULL;
	}
	//texture used to store features
	if(	_featureTex)
	{
		delete [] _featureTex;
		_featureTex =	NULL;
	}
	//texture used for multi-orientation 
	if(_orientationTex)
	{
		delete [] _orientationTex;
		_orientationTex = NULL;
	}
	int no = _octave_num* param._dog_level_num;

	//two sets of vbos used to display the features
	if(_featureDisplayVBO)
	{
		glDeleteBuffers(no, _featureDisplayVBO);
		delete [] _featureDisplayVBO;
		_featureDisplayVBO = NULL;
	}
	if( _featurePointVBO)
	{
		glDeleteBuffers(no, _featurePointVBO);
		delete [] _featurePointVBO;
		_featurePointVBO = NULL;
	}

}

void PyramidGL::DestroySharedData()
{
	//histogram reduction
	if(_histoPyramidTex)
	{
		delete[]	_histoPyramidTex;
		_hpLevelNum = 0;
		_histoPyramidTex = NULL;
	}

	//descriptor storage shared by all levels
	if(_descriptorTex)
	{
		delete [] _descriptorTex;
		_descriptorTex = NULL;
	}
	//cpu reduction buffer.
	if(_histo_buffer)
	{
		delete[] _histo_buffer;
		_histo_buffer = 0;
	}
}

void PyramidNaive::FitHistogramPyramid()
{
	GLTexImage * tex, *htex;
	int hist_level_num = _hpLevelNum - _pyramid_octave_first; 

	tex = GetBaseLevel(_octave_min , DATA_KEYPOINT) + 2;
	htex = _histoPyramidTex + hist_level_num - 1;
	int w = tex->GetImgWidth() >> 1;
	int h = tex->GetImgHeight() >> 1;

	for(int k = 0; k <hist_level_num -1; k++, htex--)
	{
		if(htex->GetImgHeight()!= h || htex->GetImgWidth() != w)
		{	
			htex->SetImageSize(w, h);
			htex->ZeroHistoMargin();
		}

		w = (w + 1)>>1; h = (h + 1) >> 1;
	}
}

void PyramidNaive::FitPyramid(int w, int h)
{
	//(w, h) <= (_pyramid_width, _pyramid_height);

	_pyramid_octave_first = 0;
	//
	_octave_num  = GlobalUtil::_octave_num_default;

	int _octave_num_max = max(1, (int) floor (log ( double(min(w, h)))/log(2.0))  -3 );

	if(_octave_num < 1 || _octave_num > _octave_num_max) 
	{
		_octave_num = _octave_num_max;
	}


	int pw = _pyramid_width>>1, ph = _pyramid_height>>1;
	while(_pyramid_octave_first + _octave_num < _pyramid_octave_num &&  
		pw >= w && ph >= h)
	{
		_pyramid_octave_first++;
		pw >>= 1;
		ph >>= 1;
	}

	for(int i = 0; i < _octave_num; i++)
	{
		GLTexImage * tex = GetBaseLevel(i + _octave_min);
		GLTexImage * aux = GetBaseLevel(i + _octave_min, DATA_KEYPOINT);
		for(int j = param._level_min; j <= param._level_max; j++, tex++, aux++)
		{
			tex->SetImageSize(w, h);
			aux->SetImageSize(w, h);
		}
		w>>=1;
		h>>=1;
	}
}
void PyramidNaive::InitPyramid(int w, int h, int ds)
{
	int wp, hp, toobig = 0;
	if(ds == 0)
	{
		_down_sample_factor = 0;
		if(GlobalUtil::_octave_min_default>=0)
		{
			wp = w >> GlobalUtil::_octave_min_default;
			hp = h >> GlobalUtil::_octave_min_default;
		}else 
		{
			wp = w << (-GlobalUtil::_octave_min_default);
			hp = h << (-GlobalUtil::_octave_min_default);
		}
		_octave_min = _octave_min_default;
	}else
	{
		//must use 0 as _octave_min; 
		_octave_min = 0;
		_down_sample_factor = ds;
		w >>= ds;
		h >>= ds;
		wp = w;
		hp = h; 

	}

	while(wp > GlobalUtil::_texMaxDim || hp > GlobalUtil::_texMaxDim)
	{
		_octave_min ++;
		wp >>= 1;
		hp >>= 1;
		toobig = 1;
	}

	if(toobig && GlobalUtil::_verbose)
	{

		std::cout<< "**************************************************************\n"
					"Image larger than allowed dimension, data will be downsampled!\n"
					"use -maxd to change the settings\n"
					"***************************************************************\n";
	}

	if( wp == _pyramid_width && hp == _pyramid_height && _allocated )
	{
		FitPyramid(wp, hp);
	}else if(GlobalUtil::_ForceTightPyramid || _allocated ==0)
	{
		ResizePyramid(wp, hp);
	}
	else if( wp > _pyramid_width || hp > _pyramid_height )
	{
		ResizePyramid(max(wp, _pyramid_width), max(hp, _pyramid_height));
		if(wp < _pyramid_width || hp < _pyramid_height)  FitPyramid(wp, hp);
	}
	else
	{
		//try use the pyramid allocated for large image on small input images
		FitPyramid(wp, hp);
	}

	//select the initial smoothing filter according to the new _octave_min
	ShaderMan::SelectInitialSmoothingFilter(_octave_min, param);
}

void PyramidNaive::ResizePyramid( int w,  int h)
{
	//
	unsigned int totalkb = 0;
	int _octave_num_new, input_sz;
	int i, j;
	GLTexImage * tex, *aux;
	//

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

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

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

	
	//compute # of octaves

	input_sz = min(w,h) ;


	_pyramid_width =  w;
	_pyramid_height =  h;

	//reset to preset parameters
	_octave_num_new  = GlobalUtil::_octave_num_default;

	if(_octave_num_new < 1) 
	{
		_octave_num_new = (int) floor (log ( double(input_sz))/log(2.0)) -3 ;
		if(_octave_num_new<1 ) _octave_num_new = 1;
	}

	if(_pyramid_octave_num != _octave_num_new)
	{
		//destroy the original pyramid if the # of octave changes
		if(_octave_num >0)
		{
			DestroyPerLevelData();
			DestroyPyramidData();
		}
		_pyramid_octave_num = _octave_num_new;
	}

	_octave_num = _pyramid_octave_num;

	int noct = _octave_num;
	int nlev = param._level_num;

	//	//initialize the pyramid
	if(_texPyramid==NULL)	_texPyramid = new GLTexImage[ noct* nlev ];
	if(_auxPyramid==NULL)	_auxPyramid = new GLTexImage[ noct* nlev ];


	tex = GetBaseLevel(_octave_min, DATA_GAUSSIAN);
	aux = GetBaseLevel(_octave_min, DATA_KEYPOINT);
	for(i = 0; i< noct; i++)
	{
		totalkb += (nlev * w * h * 16 / 1024);
		for( j = 0; j< nlev; j++, tex++)
		{
			tex->InitTexture(w, h);
			//tex->AttachToFBO(0);
		}
		//several auxilary textures are not actually required
		totalkb += ((nlev - 3) * w * h * 16 /1024);
		for( j = 0; j< nlev ; j++, aux++)
		{
			if(j < 2) continue;
			if(j >= nlev - 1) continue;
			aux->InitTexture(w, h, 0);
			//aux->AttachToFBO(0);
		}

		w>>=1;
		h>>=1;
	}

	totalkb += ResizeFeatureStorage();


	//
	_allocated = 1;

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

}


int PyramidGL::ResizeFeatureStorage()
{
	int totalkb = 0;
	if(_levelFeatureNum==NULL)	_levelFeatureNum = new int[_octave_num * param._dog_level_num];
	std::fill(_levelFeatureNum, _levelFeatureNum+_octave_num * param._dog_level_num, 0); 

	int wmax = GetBaseLevel(_octave_min)->GetDrawWidth();
	int hmax = GetBaseLevel(_octave_min)->GetDrawHeight();
	int w ,h, i;

	//use a fbo to initialize textures..
	FrameBufferObject fbo;
	
	//
	if(_histo_buffer == NULL) _histo_buffer = new float[1 << (2 + 2 * GlobalUtil::_ListGenSkipGPU)];
	//histogram for feature detection

	int num = (int)ceil(log(double(max(wmax, hmax)))/log(2.0));

	if( _hpLevelNum != num)
	{
		_hpLevelNum = num;
		if(GlobalUtil::_ListGenGPU)
		{
			if(_histoPyramidTex ) delete [] _histoPyramidTex;
			_histoPyramidTex = new GLTexImage[_hpLevelNum];
			w = h = 1 ;
			for(i = 0; i < _hpLevelNum; i++)
			{
				_histoPyramidTex[i].InitTexture(w, h, 0);
				_histoPyramidTex[i].AttachToFBO(0);
				w<<=1;
				h<<=1;
			}
		}
	}

	// (4 ^ (_hpLevelNum) -1 / 3) pixels
	if(GlobalUtil::_ListGenGPU) totalkb += (((1 << (2 * _hpLevelNum)) -1) / 3 * 16 / 1024);



	//initialize the feature texture

	int idx = 0, n = _octave_num * param._dog_level_num;
	if(_featureTex==NULL)	_featureTex = new GLTexImage[n];
	if(GlobalUtil::_MaxOrientation >1 && GlobalUtil::_OrientationPack2==0)
	{
		if(_orientationTex== NULL)		_orientationTex = new GLTexImage[n];
	}


	for(i = 0; i < _octave_num; i++)
	{
		GLTexImage * tex = GetBaseLevel(i+_octave_min);
		int fmax = int(tex->GetImgWidth()*tex->GetImgHeight()*GlobalUtil::_MaxFeaturePercent);
		int fw, fh;
		//
		if(fmax > GlobalUtil::_MaxLevelFeatureNum) fmax = GlobalUtil::_MaxLevelFeatureNum;
		else if(fmax < 32) fmax = 32;	//give it at least a space of 32 feature

		GetTextureStorageSize(fmax, fw, fh);
		
		for(int j = 0; j < param._dog_level_num; j++, idx++)
		{