siftgpu.cpp

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

		//Loop between verbose level 0, 1, 2
		if(GlobalUtil::_verbose)
		{
			GlobalUtil::_verbose  = GlobalUtil::_timingS;
			GlobalUtil::_timingS = 0;
			if(GlobalUtil::_verbose ==0 && GlobalUtil::_UseSiftGPUEX)
				std::cout << "Console ouput disabled, press Q/V to enable\n\n";
		}else
		{
			GlobalUtil::_verbose = 1;
			GlobalUtil::_timingS = 1;
		}
	}else if(verbose == -2)
	{
		//trick for disabling all output (still keeps the timing level)
		GlobalUtil::_verbose = 0;
		GlobalUtil::_timingS = 1;
	}else 
	{
		GlobalUtil::_verbose = verbose>0;
		GlobalUtil::_timingS = verbose>1;
	}

}


SiftParam::SiftParam()
{

	_level_min = -1;
	_dog_level_num  = 3;
	_level_max = 0;
	_sigma0 = 0;
	_sigman = 0;
	_edge_threshold = 0;
	_dog_threshold =  0;


}

float SiftParam::GetInitialSmoothSigma(int octave_min)
{
	float	sa = _sigma0 * powf(2.0f, float(_level_min)/float(_dog_level_num)) ; 
	float   sb = _sigman / powf(2.0f,  float(octave_min)) ;//
	float sigma_skip0 = sa>sb+ 0.001?sqrt(sa*sa - sb*sb): 0.0f;
	return sigma_skip0; 
}

void SiftParam::ParseSiftParam()
{ 

	if(_dog_level_num ==0) _dog_level_num = 3;
	if(_level_max ==0) _level_max = _dog_level_num + 1;
	if(_sigma0 ==0.0f) _sigma0 = 1.6f * powf(2.0f, 1.0f / _dog_level_num) ;
	if(_sigman == 0.0f) _sigman = 0.5f;


	_level_num = _level_max -_level_min + 1;

	_level_ds  = _level_min + _dog_level_num;
	if(_level_ds > _level_max ) _level_ds = _level_max ;


	///
	float _sigmak = powf(2.0f, 1.0f / _dog_level_num) ;
	float dsigma0 = _sigma0 * sqrt (1.0f - 1.0f / (_sigmak*_sigmak) ) ;
	float sa, sb;

 
	sa = _sigma0 * powf(_sigmak, (float)_level_min) ; 
	sb = _sigman / powf(2.0f,   (float)GlobalUtil::_octave_min_default) ;//

	_sigma_skip0 = sa>sb+ 0.001?sqrt(sa*sa - sb*sb): 0.0f;

    sa = _sigma0 * powf(_sigmak, float(_level_min )) ;
    sb = _sigma0 * powf(_sigmak, float(_level_ds - _dog_level_num)) ;

	_sigma_skip1 = sa>sb + 0.001? sqrt(sa*sa - sb*sb): 0.0f;

	_sigma_num = _level_max - _level_min;
	_sigma = new float[_sigma_num];

	for(int i = _level_min + 1; i <= _level_max; i++)
	{
		_sigma[i-_level_min -1] =  dsigma0 * powf(_sigmak, float(i)) ;
	}

	if(_dog_threshold ==0)	_dog_threshold      = 0.02f / _dog_level_num ;
	if(_edge_threshold==0) _edge_threshold		= 10.0f;
}


void SiftGPUEX::DisplayOctave(void (*UseDisplayShader)(), int i)
{
	if(_pyramid == NULL)return;
	const int grid_sz = (int)ceil(_level_num/2.0);
	double scale = 1.0/grid_sz ;
	int gx=0, gy=0, dx, dy;

	if(_pyramid->_octave_min >0) scale *= (1<<_pyramid->_octave_min);
	else if(_pyramid->_octave_min < 0) scale /= (1<<(-_pyramid->_octave_min));


	i = i% _pyramid->_octave_num;  //
	if(i<0 ) i+= _pyramid->_octave_num;

	scale *= ( 1<<(i));




	UseDisplayShader();

	glPushMatrix();
	glScaled(scale, scale, scale);
	for(int level = _level_min; level<= _level_max; level++)
	{
		GLTexImage * tex = _pyramid->GetLevelTexture(i+_pyramid->_octave_min, level);

		dx = tex->GetImgWidth();
		dy = tex->GetImgHeight();

		glPushMatrix();

		glTranslated(dx*gx, dy*gy, 0);

		tex->BindTex();

		tex->DrawImage();
		tex->UnbindTex();

		glPopMatrix();

		gx++;
		if(gx>=grid_sz) 
		{
			gx =0;
			gy++;
		}

	}

	glPopMatrix();
	ShaderMan::UnloadProgram();
}

void SiftGPUEX::DisplayPyramid( void (*UseDisplayShader)(), int dataName, int nskip1, int nskip2)
{

	if(_pyramid == NULL)return;
	int grid_sz = (_level_num -nskip1 - nskip2);
	if(grid_sz > 4) grid_sz = (int)ceil(grid_sz*0.5);
	double scale = 1.0/grid_sz;
	int stepx = 0, stepy = 0, dx, dy=0, nstep;

	if(_pyramid->_octave_min >0) scale *= (1<<_pyramid->_octave_min);
	else if(_pyramid->_octave_min < 0) scale /= (1<<(-_pyramid->_octave_min));


	glPushMatrix();
	glScaled(scale, scale, scale);

	for(int i = _pyramid->_octave_min; i < _pyramid->_octave_min+_pyramid->_octave_num; i++)
	{
	
		nstep = i==_pyramid->_octave_min? grid_sz: _level_num;
		dx = 0;
		UseDisplayShader();
		for(int j = _level_min + nskip1; j <= _level_max-nskip2; j++)
		{
			GLTexImage * tex = _pyramid->GetLevelTexture(i, j, dataName);
			if(tex->GetImgWidth() == 0 || tex->GetImgHeight() == 0) continue;
			stepx = tex->GetImgWidth();
			stepy = tex->GetImgHeight();	
			////
			if(j == _level_min + nskip1 + nstep)
			{
				dy += stepy;
				dx = 0;
			}

			glPushMatrix();
			glTranslated(dx, dy, 0);
			tex->BindTex();
			tex->DrawImage();
			tex->UnbindTex();
			glPopMatrix();

			dx += stepx;

		}

		ShaderMan::UnloadProgram();

		dy+= stepy;
	}

	glPopMatrix();
}


void SiftGPUEX::DisplayLevel(void (*UseDisplayShader)(), int i)
{
	if(_pyramid == NULL)return;

	i = i%(_level_num * _pyramid->_octave_num);
	if (i<0 ) i+= (_level_num * _pyramid->_octave_num);
	int octave = _pyramid->_octave_min + i/_level_num;
	int level  = _level_min + i%_level_num;
	double scale = 1.0;

	if(octave >0) scale *= (1<<octave);
	else if(octave < 0) scale /= (1<<(-octave));

	GLTexImage * tex = _pyramid->GetLevelTexture(octave, level);

	UseDisplayShader();

	glPushMatrix();
	glScaled(scale, scale, scale);
	tex->BindTex();
	tex->DrawImage();
	tex->UnbindTex();
	glPopMatrix();
	ShaderMan::UnloadProgram();
}

void SiftGPUEX::DisplaySIFT()
{
	if(_pyramid == NULL) return;
	if(_view_debug)
	{
		DisplayDebug();
		return;
	}
	switch(_view)
	{
	case 0:
		DisplayInput();
		DisplayFeatureBox(_sub_view);
		break;
	case 1:
		DisplayPyramid(ShaderMan::UseShaderDisplayGaussian, SiftPyramid::DATA_GAUSSIAN);
		break;
	case 2:
		DisplayOctave(ShaderMan::UseShaderDisplayGaussian, _sub_view);	
		break;
	case 3:
		DisplayLevel(ShaderMan::UseShaderDisplayGaussian, _sub_view);
		break;
	case 4:
		DisplayPyramid(ShaderMan::UseShaderDisplayDOG, SiftPyramid::DATA_DOG, 1);
		break;
	case 5:
		DisplayPyramid(ShaderMan::UseShaderDisplayGrad, SiftPyramid::DATA_GRAD, 1);
		break;
	case 6:
		DisplayPyramid(ShaderMan::UseShaderDisplayDOG, SiftPyramid::DATA_DOG,2, 1);
		DisplayPyramid(ShaderMan::UseShaderDisplayKeypoints, SiftPyramid::DATA_KEYPOINT, 2,1);
	}
}


void SiftGPUEX::SetView(int view, int sub_view, char *title)
{
	const char* view_titles[] =
	{
		"Original Image",
		"Gaussian Pyramid",
		"Octave Images",
		"Level Image",
		"Difference of Gaussian",
		"Gradient",
		"Keypoints"
	};
	const int view_num = 7;
	_view = view % view_num;
	if(_view <0) _view +=view_num;
	_sub_view = sub_view;

	if(_view_debug)
		strcpy(title, "Debug...");
	else
		strcpy(title, view_titles[_view]);

}


void SiftGPU::PrintUsage()
{
	std::cout
	<<"SiftGPU Usage:\n"
	<<"-h -help       : Help\n"
	<<"-i <strings>   : Filename(s) of the input image(s)\n"
	<<"-il <string>   : Filename of an image list file\n"
	<<"-o <string>    : Where to save SIFT features\n"
	<<"-f <float>     : Filter width factor; Width will be 2*factor+1 (default : 4.0)\n"
	<<"-w  <float>    : Orientation sample window factor (default: 2.0)\n"
	<<"-dw <float>  * : Descriptor grid size factor (default : 3.0)\n"
	<<"-fo <int>    * : First octave to detect DOG keypoints(default : 0)\n"
	<<"-no <int>      : Maximum number of Octaves (default : no limit)\n"
	<<"-d <int>       : Number of DOG levels in an octave (default : 3)\n"
	<<"-t <float>     : DOG threshold (default : 0.02/3)\n"
	<<"-e <float>     : Edge Threshold (default : 10.0)\n"
	<<"-m  <int=2>    : Multi Feature Orientations (default : 1)\n"
	<<"-m2p           : 2 Orientations packed as one float\n"
	<<"-s  <int=1>    : Sub-Pixel, Sub-Scale Localization, Multi-Refinement(num)\n"
	<<"-lcpu -lc <int>: CPU/GPU mixed Feature List Generation (defaut : 6)\n"
	<<"                 Use GPU first, and use CPU when reduction size <= pow(2,num)\n"
	<<"                 When <num> is missing or equals -1, no GPU will be used\n"
	<<"-noprep        : Upload raw data to GPU (default: RGB->LUM and down-sample on CPU)\n"
	<<"-sd            : Skip descriptor computation if specified\n"
	<<"-unn    *      : Write unnormalized descriptor if specified\n"
	<<"-b      *      : Write binary sift file if specified\n"
	<<"-fs <int>      : Block Size for freature storage <default : 4>\n"
	<<"-glsl          : Use GLSL SIFTGPU instead of CG (default : CG)\n"
	<<"-tight         : Automatically resize pyramid to fit new images tightly\n"
	<<"-p  WxH        : Inititialize the pyramids to contain image of WxH (eg -p 1024x768)\n"
	<<"-lm  <int>     : Maximum feature count for a level (for pre-allocation)\n"
	<<"-lmp <float>   : Maximum percent of pixels as features (for pre-allocaton)\n"
	<<"-v <int>       : Level of timing details. Same as calling Setverbose() function\n"
	<<"-loweo         : (0,0) at center of top-left pixel (defaut: corner)\n"
	<<"-maxd <int> *  : Max working dimension (default : 2560 (unpacked) / 3200 (packed))\n"
	<<"-exit          : Exit program after processing the input image\n"
	<<"-unpack        : Use the old unpacked implementation\n"
	<<"-di            : Use dynamic array indexing if available (defualt : no)\n"
	<<"                 It could make computation faster on cards like GTX 280\n"
	<<"-fastmath      : specify -fastmath to cg compiler (Not much difference.)\n"
	<<"-ofix     *    : use 0 as feature orientations.\n"
	<<"-ofix-not *	  : disable -ofix.\n"
	<<"------parameters marked with * can be changed after initialization---------\n"
	<<"\n";
}
void SiftGPU::ParseParam(int argc, char **argv)
{
	char* arg, *param;
	char* opt;
	int  HelpPrinted = 0, setMaxD = 0;
	int  i = 0;
	for( i = 0; i< argc; i++)
	{
		arg = argv[i];
		if(arg[0]!='-')continue;
		opt = arg+1;
		param = argv[i+1];
		if(_stricmp(opt, "h")==0 || _stricmp(opt,"help")==0)
		{
			HelpPrinted = 1;
			PrintUsage();
		}else if(_stricmp(opt, "glsl")==0)
		{
			GlobalUtil::_UseGLSL = 1;
		}else if(_stricmp(opt, "cg")==0)
		{
#if !defined(SIFTGPU_NO_CG)
			GlobalUtil::_UseGLSL = 0;
#endif
		}else if(_stricmp(opt, "cuda")==0)
		{
			GlobalUtil::_UseCUDA = 1;
			GlobalUtil::_usePackedTex = 0;
		}else if(_stricmp(opt, "pack")==0)
		{
			GlobalUtil::_usePackedTex = 1;
		}else if(_stricmp(opt, "unpack")==0)
		{
			GlobalUtil::_usePackedTex = 0;
		}else if(_stricmp(opt, "lcpu")==0||_stricmp(opt, "lc")==0)
		{
			int gskip = -1;
			if(i+1 <argc)	sscanf(param, "%d", &gskip);
			if(gskip >= 0)
			{
				GlobalUtil::_ListGenSkipGPU = gskip;
			}else
			{
				GlobalUtil::_ListGenGPU = 0;
			}
		}else if(_stricmp(opt, "prep")==0)
		{
			GlobalUtil::_PreProcessOnCPU = 1;
		}else if(_stricmp(opt, "noprep")==0)
		{
			GlobalUtil::_PreProcessOnCPU = 0;
		}else  if(_stricmp(opt, "fbo1")==0)
		{
			FrameBufferObject::UseSingleFBO =1;

		}else  if(_stricmp(opt, "fbos")==0)
		{
			FrameBufferObject::UseSingleFBO = 0;
		}
		else if(_stricmp(opt, "sd")==0)
		{
			GlobalUtil::_DescriptorPPT =0;
		}else if(_stricmp(opt, "unn")==0)
		{
			GlobalUtil::_NormalizedSIFT =0;
		}else if(_stricmp(opt, "b")==0)
		{
			GlobalUtil::_BinarySIFT = 1;
		}else if(_stricmp(opt, "tight")==0)
		{
			GlobalUtil::_ForceTightPyramid = 1;
		}else if(_stricmp(opt, "exit")==0)
		{
			GlobalUtil::_ExitAfterSIFT = 1;
		}else if(_stricmp(opt, "di")==0)
		{
			GlobalUtil::_UseDynamicIndexing = 1;
		}else if(_stricmp(opt, "sign")==0)
		{
			GlobalUtil::_KeepExtremumSign = 1;
		}else if(_stricmp(opt, "ov292")==0)
		{
			//for compatibility with old version//
			GlobalUtil::_GradientLevelOffset = 2; 
			GlobalUtil::_OrientationGaussianFactor = 4.5;
			GlobalUtil::_OrientationWindowFactor = 1.0;

		}else if(_stricmp(opt, "m")==0 || _stricmp(opt, "mo")==0)
		{
			int mo = 2; //default multi-orientation
			if(i+1 <argc)	sscanf(param, "%d", &mo);
			//at least two orientation
			GlobalUtil::_MaxOrientation = min(max(1, mo), 4);
		}else if(_stricmp(opt, "m2p") == 0)
		{
			GlobalUtil::_MaxOrientation = 2;
			GlobalUtil::_OrientationPack2 = 1;
		}else if(_stricmp(opt, "s") ==0)
		{
			int sp = 1; //default refinement
			if(i+1 <argc)	sscanf(param, "%d", &sp);
			//at least two orientation