siftmatch.cpp

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

			"		buf = texture2DRect(tex2, tpos.yw);\n"
			"		data1 = texture2DRect(tex1, tpos.xz);\n"
			"		val += (data1 * buf);\n"
			"		tpos.xy = tpos.xy + vec2(1.0, 1.0);\n"
			"	}\n"
			"	const float factor = 0.248050689697265625; \n"
			"	gl_FragColor =vec4(dot(val, vec4(factor)), index,  0);\n"
			"}"
		<<	'\0';

	s_multiply = program= new ProgramGLSL(buffer); 

	_param_multiply_tex1 = glGetUniformLocation(*program, "tex1");
	_param_multiply_tex2 = glGetUniformLocation(*program, "tex2");
	_param_multiply_size = glGetUniformLocation(*program, "size");

	out.seekp(ios::beg);
	out <<  "#pragma optionNV(ifcvt none)\n"
			"#pragma optionNV(unroll all)\n"
	        "#define SIFT_PER_STRIPE " << _sift_per_stripe << ".0\n" 
			"#define PIXEL_PER_SIFT " << _pixel_per_sift << "\n"
			"uniform sampler2DRect tex1, tex2;\n"
			"uniform sampler2DRect texL1;\n"
			"uniform sampler2DRect texL2; \n"
			"uniform mat3 H; \n"
			"uniform mat3 F; \n"
			"uniform vec4	size; \n"
			"void main()		\n"
		    "{\n"
		<<	"   vec4 val = vec4(0, 0, 0, 0), data1, buf;\n"
			"   vec2 index = gl_FragCoord.yx; \n"
			"   vec2 stripe_size = size.xy * SIFT_PER_STRIPE;\n"
			"	vec2 temp_div1 = index / stripe_size;\n"
			"   vec2 stripe_index = floor(temp_div1);\n"
			"   index = floor(stripe_size * (temp_div1 - stripe_index));\n"
			"	vec2 temp_div2 = index  * vec2(1.0/ float(SIFT_PER_STRIPE));\n"
			"	vec2 temp_floor2 = floor(temp_div2);\n"
			"   vec2 index_v = temp_floor2 + vec2(0.5);\n "
			"   vec2 index_h = vec2(SIFT_PER_STRIPE)* (temp_div2 - temp_floor2);\n"
			
			//read feature location data
			"   vec4 tlpos = vec4((index_h + stripe_index * vec2(SIFT_PER_STRIPE)) + 0.5, index_v);\n"
			"   vec3 loc1 = vec3(texture2DRect(texL1, tlpos.xz).xw, 1.0);\n"
			"   vec3 loc2 = vec3(texture2DRect(texL2, tlpos.yw).xw, 1.0);\n"
			
			//check the guiding homography
			"   vec3 hxloc1 = H* loc1;\n"
			"   vec2 diff = abs(loc2.xy- (hxloc1.xy/hxloc1.z));\n"
			"   float disth = max(diff.x, diff.y);\n"
			"   if(disth > size.z ) {gl_FragColor = vec4(0, index, 0); return;}\n"

			//check the guiding fundamental 
			"   vec3 fx1 = (F * loc1), ftx2 = (loc2 * F);\n"
			"   float x2tfx1 = dot(loc2, fx1);\n"
			"   vec4 temp = vec4(fx1.xy, ftx2.xy); \n"
			"   float sampson_error = (x2tfx1 * x2tfx1) / dot(temp, temp);\n"
			"   if(sampson_error > size.w) {gl_FragColor = vec4(0, index, 0); return;}\n"

			//compare feature descriptor
			"   vec2 tx = (index_h + stripe_index * SIFT_PER_STRIPE)* vec2(PIXEL_PER_SIFT) + 0.5;\n"
			"   vec2 tpos1, tpos2; \n"
			"	vec4 tpos = vec4(tx, index_v);\n"
			"   for(int i = 0; i < PIXEL_PER_SIFT; ++i){\n"
			"		buf = texture2DRect(tex2, tpos.yw);\n"
			"		data1 = texture2DRect(tex1, tpos.xz);\n"
			"		val += data1 * buf;\n"
			"		tpos.xy = tpos.xy + vec2(1.0, 1.0);\n"
			"	}\n"
			"	const float factor = 0.248050689697265625; \n"
			"	gl_FragColor =vec4(dot(val, vec4(factor)), index,  0);\n"
			"}"
		<<	'\0';

	s_guided_mult = program= new ProgramGLSL(buffer);

	_param_guided_mult_tex1 = glGetUniformLocation(*program, "tex1");
	_param_guided_mult_tex2= glGetUniformLocation(*program, "tex2");
	_param_guided_mult_texl1 = glGetUniformLocation(*program, "texL1");
	_param_guided_mult_texl2 = glGetUniformLocation(*program, "texL2");
	_param_guided_mult_h = glGetUniformLocation(*program, "H");
	_param_guided_mult_f = glGetUniformLocation(*program, "F");
	_param_guided_mult_param = glGetUniformLocation(*program, "size");

	//row max
	out.seekp(ios::beg);
	out <<	"#define BLOCK_WIDTH 16.0\n"
			"uniform sampler2DRect tex;	uniform vec3 param;\n"
			"void main ()\n"
			"{\n"
			"	float index = gl_FragCoord.x + floor(gl_FragCoord.y) * BLOCK_WIDTH; \n"
			"	vec2 bestv = vec2(-1.0); float imax = -1.0;\n"
			"	for(float i = 0.0; i < param.x; i ++){\n "
			"		float v = texture2DRect(tex, vec2(i + 0.5, index)).r; \n"
			"		imax = v > bestv.r ? i : imax; \n "
			"		bestv  = v > bestv.r? vec2(v, bestv.r) : max(bestv, vec2(v));\n "
			"	}\n"
			"	bestv = acos(min(bestv, 1.0));\n"
			"	if(bestv.x >= param.y || bestv.x >= param.z * bestv.y) imax = -1.0;\n"
			"	gl_FragColor = vec4(imax, bestv, index);\n"
			"}"
		<<  '\0';
	s_row_max = program= new ProgramGLSL(buffer); 
	_param_rowmax_param = glGetUniformLocation(*program, "param");

	out.seekp(ios::beg);
	out <<	"#define BLOCK_WIDTH 16.0\n"
			"uniform sampler2DRect tex; uniform vec3 param;\n"
			"void main ()\n"
			"{\n"
			"	float index = gl_FragCoord.x + floor(gl_FragCoord.y) * BLOCK_WIDTH; \n"
			"	vec2 bestv = vec2(-1.0); float imax = -1.0;\n"
			"	for(float i = 0.0; i < param.x; i ++){\n "
			"		float v = texture2DRect(tex, vec2(index, i + 0.5)).r; \n"
			"		imax = (v > bestv.r)? i : imax; \n "
			"		bestv  = v > bestv.r? vec2(v, bestv.r) : max(bestv, vec2(v));\n "
			"	}\n"
			"	bestv = acos(min(bestv, 1.0));\n"
			"	if(bestv.x >= param.y || bestv.x >= param.z * bestv.y) imax = -1.0;\n"
			"	gl_FragColor = vec4(imax, bestv, index);\n"
			"}"
		<<  '\0';
	s_col_max = program =new ProgramGLSL(buffer); 
	_param_colmax_param = glGetUniformLocation(*program, "param");


}

int  SiftMatchGL::GetGuidedSiftMatch(int max_match, int match_buffer[][2], float H[3][3], float F[3][3],
									 float distmax, float ratiomax, float hdistmax, float fdistmax, int mbm)
{

	int dw = _num_sift[1];
	int dh = _num_sift[0]; 
	if(_initialized ==0) return 0;
	if(dw <= 0 || dh <=0) return 0;
	if(_have_loc[0] == 0 || _have_loc[1] == 0) return 0;

	FrameBufferObject fbo;
	glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
	_texDot.SetImageSize(dw, dh);


	//data
	_texDot.AttachToFBO(0);
	_texDot.FitTexViewPort();
	glActiveTexture(GL_TEXTURE0);
	_texDes[0].BindTex();
	glActiveTexture(GL_TEXTURE1);
	_texDes[1].BindTex();
	glActiveTexture(GL_TEXTURE2);
	_texLoc[0].BindTex();
	glActiveTexture(GL_TEXTURE3);
	_texLoc[1].BindTex();

	//multiply the descriptor matrices
	s_guided_mult->UseProgram();


	//set parameters glsl/cg
	float dot_param[4] = {(float)_texDes[0].GetDrawHeight(), (float) _texDes[1].GetDrawHeight(), hdistmax, fdistmax};
#if !defined(SIFTGPU_NO_CG)
	if(_use_glsl)
	{
#endif
		glUniform1i(_param_guided_mult_tex1, 0);
		glUniform1i(_param_guided_mult_tex2, 1);
		glUniform1i(_param_guided_mult_texl1, 2);
		glUniform1i(_param_guided_mult_texl2, 3);
		glUniformMatrix3fv(_param_guided_mult_h, 1, GL_TRUE, H[0]);
		glUniformMatrix3fv(_param_guided_mult_f, 1, GL_TRUE, F[0]);
		glUniform4fv(_param_guided_mult_param, 1, dot_param);
#if !defined(SIFTGPU_NO_CG)
	}else
	{

		cgGLSetTextureParameter(_param_guided_mult_tex1, _texDes[0]);	
		cgGLEnableTextureParameter(_param_guided_mult_tex1);
		cgGLSetTextureParameter(_param_guided_mult_tex2, _texDes[1]);	
		cgGLEnableTextureParameter(_param_guided_mult_tex2);
		cgGLSetTextureParameter(_param_guided_mult_texl1, _texLoc[0]);	
		cgGLEnableTextureParameter(_param_guided_mult_texl1);
		cgGLSetTextureParameter(_param_guided_mult_texl2, _texLoc[1]);	
		cgGLEnableTextureParameter(_param_guided_mult_texl2);
		cgSetMatrixParameterfr(_param_guided_mult_h, H[0]);
		cgSetMatrixParameterfr(_param_guided_mult_f, F[0]);
		cgGLSetParameter4fv(_param_guided_mult_param, dot_param);
	}
#endif


	_texDot.DrawQuad();

	GLTexImage::UnbindMultiTex(4);

	return GetBestMatch(max_match, match_buffer, distmax, ratiomax, mbm);
}

int SiftMatchGL::GetBestMatch(int max_match, int match_buffer[][2], float distmax, float ratiomax, int mbm)
{

	glActiveTexture(GL_TEXTURE0);
	_texDot.BindTex();

	//readback buffer
	sift_buffer.resize(_num_sift[0] + _num_sift[1] + 16);
	float * buffer1 = &sift_buffer[0], * buffer2 = &sift_buffer[_num_sift[0]];

	//row max
	_texMatch[0].AttachToFBO(0);
	_texMatch[0].SetImageSize(16, ( _num_sift[0] + 15) / 16);
	_texMatch[0].FitTexViewPort();

	///set parameter glsl/cg
	s_row_max->UseProgram();
#if !defined(SIFTGPU_NO_CG)
	if(_use_glsl) glUniform3f(_param_rowmax_param, (float)_num_sift[1], distmax, ratiomax);
	else cgGLSetParameter3f(_param_rowmax_param, (float)_num_sift[1], distmax, ratiomax);
#else
	glUniform3f(_param_rowmax_param, (float)_num_sift[1], distmax, ratiomax);
#endif

	_texMatch[0].DrawQuad();
	glReadPixels(0, 0, 16, (_num_sift[0] + 15)/16, GL_RED, GL_FLOAT, buffer1);

	//col max
	if(mbm)
	{
		_texMatch[1].AttachToFBO(0);
		_texMatch[1].SetImageSize(16, (_num_sift[1] + 15) / 16);
		_texMatch[1].FitTexViewPort();
		//set parameter glsl/cg
		s_col_max->UseProgram();
#if !defined(SIFTGPU_NO_CG)
		if(_use_glsl) glUniform3f(_param_rowmax_param, (float)_num_sift[0], distmax, ratiomax);
		else cgGLSetParameter3f(_param_colmax_param, (float)_num_sift[0], distmax, ratiomax);
#else
		glUniform3f(_param_rowmax_param, (float)_num_sift[0], distmax, ratiomax);
#endif
		_texMatch[1].DrawQuad();
		glReadPixels(0, 0, 16, (_num_sift[1] + 15) / 16, GL_RED, GL_FLOAT, buffer2);
	}


	//unload
#if !defined(SIFTGPU_NO_CG)
	if(_use_glsl)
	{
#endif
		glUseProgram(0);
#if !defined(SIFTGPU_NO_CG)
	}else
	{
		cgGLUnbindProgram(ProgramCG::_FProfile);
		cgGLDisableProfile(ProgramCG::_FProfile);
	}
#endif
	GLTexImage::UnbindMultiTex(2);
	GlobalUtil::CleanupOpenGL();

	//write back the matches
	int nmatch = 0, j ;
	for(int i = 0; i < _num_sift[0] && nmatch < max_match; ++i)
	{
		j = int(buffer1[i]);
		if( j>= 0 && (!mbm ||int(buffer2[j]) == i))
		{
			match_buffer[nmatch][0] = i;
			match_buffer[nmatch][1] = j;
			nmatch++;
		}
	}
	return nmatch;
}

int  SiftMatchGL::GetSiftMatch(int max_match, int match_buffer[][2], float distmax, float ratiomax, int mbm)
{
	int dw = _num_sift[1];
	int dh =  _num_sift[0]; 
	if(_initialized ==0) return 0;
	if(dw <= 0 || dh <=0) return 0;

	FrameBufferObject fbo;
	glDrawBuffer(GL_COLOR_ATTACHMENT0_EXT);
	_texDot.SetImageSize(dw, dh);

	//data
	_texDot.AttachToFBO(0);
	_texDot.FitTexViewPort();
	glActiveTexture(GL_TEXTURE0);
	_texDes[0].BindTex();
	glActiveTexture(GL_TEXTURE1);
	_texDes[1].BindTex();

	//////////////////	
	//multiply the descriptor matrices
	s_multiply->UseProgram();
	//set parameters
	float heights[2] = {(float)_texDes[0].GetDrawHeight(), (float)_texDes[1].GetDrawHeight()};

#if !defined(SIFTGPU_NO_CG)
	if(_use_glsl)
	{
#endif
		glUniform1i(_param_multiply_tex1, 0);
		glUniform1i(_param_multiply_tex2 , 1);
		glUniform2fv(_param_multiply_size, 1, heights);

#if !defined(SIFTGPU_NO_CG)
	}else
	{
		cgGLSetTextureParameter(_param_multiply_tex1, _texDes[0]);	
		cgGLEnableTextureParameter(_param_multiply_tex1);
		cgGLSetTextureParameter(_param_multiply_tex2, _texDes[1]);	
		cgGLEnableTextureParameter(_param_multiply_tex2);
		cgGLSetParameter2fv(_param_multiply_size, heights);
	}
#endif

	_texDot.DrawQuad();

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GlobalUtil::_texTarget, 0);

	return GetBestMatch(max_match, match_buffer, distmax, ratiomax, mbm);
}


int SiftMatchGPU::CreateContextGL()
{
	//use GLUT to create an OpenGL Context?
	if(!GlobalUtil::CreateWindowGLUT()) return 0;
	return VerifyContextGL();
}


int SiftMatchGPU::VerifyContextGL()
{
	if(__matcher) return GlobalUtil::_GoodOpenGL;
	else if(__language < 0 || __max_sift < 0)
	{
		__matcher = NULL;
		return 0;
	}else
	{
#ifdef CUDA_SIFTGPU_ENABLED
		if((__language == 0 && GlobalUtil::_UseCUDA) || __language == 3)
		{
			__language = 3;
			__matcher = new SiftMatchCU(__max_sift);
		}else
#else
		if((__language == 0 && GlobalUtil::_UseCUDA) || __language == 3)
		{
			std::cerr	<< "---------------------------------------------------------------------------\n"
						<< "CUDA is not supported in this binary! To enable CUDA implementation, please\n" 
						<< "use SiftGPU_CUDA_Enable Project for VS2005 or define CUDA_SIFTGPU_ENABLED\n"
						<< "----------------------------------------------------------------------------\n";
		}
#endif
		if((__language == 0 && GlobalUtil::_UseGLSL) || __language == 2)
		{
			__language = 2;
			__matcher = new SiftMatchGL(__max_sift, 1);
		}else
		{
			__language = 1;
			__matcher = new SiftMatchGL(__max_sift, 0);
		}
		
	}
	__matcher->InitSiftMatch();
	return GlobalUtil::_GoodOpenGL;
}

void* SiftMatchGPU::operator new (size_t  size){
  void * p = malloc(size);
  if (p == 0)  
  {
	  const std::bad_alloc ba;
	  throw ba; 
  }
  return p; 
}



SiftMatchGPU::SiftMatchGPU(int max_sift)
{
	__max_sift = max_sift;
	__language = 0;
	__matcher = NULL;
}

void SiftMatchGPU::SetLanguage(int gpu_language)
{
	if(__matcher) return;
	__language = gpu_language;
}

void SiftMatchGPU::SetMaxSift(int max_sift)
{
	if(__matcher)	__matcher->SetMaxSift(max_sift);
	else __max_sift = max_sift;
}

SiftMatchGPU::~SiftMatchGPU()
{
	if(__matcher) delete __matcher;
}

void SiftMatchGPU::SetDescriptors(int index, int num, const unsigned char* descriptors, int id)
{
	__matcher->SetDescriptors(index, num,  descriptors, id);
}

void SiftMatchGPU::SetDescriptors(int index, int num, const float* descriptors, int id)
{
	__matcher->SetDescriptors(index, num, descriptors, id);
}

void SiftMatchGPU::SetFeautreLocation(int index, const float* locations, int gap)
{
	__matcher->SetFeautreLocation(index, locations, gap);

}
int  SiftMatchGPU::GetGuidedSiftMatch(int max_match, int match_buffer[][2], float H[3][3], float F[3][3], 
				float distmax, float ratiomax, float hdistmax, float fdistmax, int mutual_best_match)
{
	if(H == NULL && F == NULL)
	{
		return __matcher->GetSiftMatch(max_match, match_buffer, distmax, ratiomax, mutual_best_match);
	}else
	{
		float Z[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, ti = (1.0e+20F);

		return __matcher->GetGuidedSiftMatch(max_match, match_buffer, H? H : Z, F? F : Z,
			distmax, ratiomax, H? hdistmax: ti,  F? fdistmax: ti, mutual_best_match);
	}
}
int  SiftMatchGPU::GetSiftMatch(int max_match, int match_buffer[][2], float distmax, float ratiomax, int mutual_best_match)
{

/*	timer.StartTimer("Sift Match 100");
	for(int i = 0;i < 100; ++i)__matcher->GetSiftMatch(max_match, match_buffer, distmax, ratiomax, mutual_best_match);
	timer.StopTimer();*/

	return __matcher->GetSiftMatch(max_match, match_buffer, distmax, ratiomax, mutual_best_match);
}

SiftMatchGPU* CreateNewSiftMatchGPU(int max_sift)
{
	return new SiftMatchGPU(max_sift);
}