siftmatch.cpp

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

////////////////////////////////////////////////////////////////////////////
//	File:		SiftMatch.cpp
//	Author:		Changchang Wu
//	Description :	implementation of SiftMatchGPU and SiftMatchGL
//
//
//	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 <strstream>
#include <algorithm>
using namespace std;
#include "GlobalUtil.h"

#if !defined(SIFTGPU_NO_CG)
#include "ProgramCG.h"
#endif

#include "ProgramGLSL.h"
#include "GLTexImage.h"
#include "SiftGPU.h"
#include "SiftMatch.h"
#include "FrameBufferObject.h"

#if (!defined(_MSC_VER) ||_MSC_VER >= 1400) && defined(CUDA_SIFTGPU_ENABLED)
#include "CuTexImage.h"
#include "SiftMatchCU.h"
#endif


SiftMatchGL::SiftMatchGL(int max_sift, int use_glsl): SiftMatchGPU()
{
	s_multiply = s_col_max = s_row_max = s_guided_mult = NULL;
	_num_sift[0] = _num_sift[1] = 0;
	_id_sift[0] = _id_sift[1] = 0;
	_have_loc[0] = _have_loc[1] = 0;
	_max_sift = max_sift <=0 ? 4096 : ((max_sift + 31)/ 32 * 32) ; 
	_pixel_per_sift = 32; //must be 32
	_sift_num_stripe = 1; 
	_sift_per_stripe = 1;
	_sift_per_row = _sift_per_stripe * _sift_num_stripe;
#if !defined(SIFTGPU_NO_CG)
	_use_glsl = use_glsl;
#else
	_use_glsl = 1;
#endif
	_initialized = 0;
}

SiftMatchGL::~SiftMatchGL()
{
	if(s_multiply) delete s_multiply;
	if(s_guided_mult) delete s_guided_mult;
	if(s_col_max) delete s_col_max;
	if(s_row_max) delete s_row_max;
}

void SiftMatchGL::SetMaxSift(int max_sift)
{
	
	max_sift = ((max_sift + 31)/32)*32;
	if(max_sift > GlobalUtil::_texMaxDimGL) max_sift = GlobalUtil::_texMaxDimGL;
	if(max_sift > _max_sift)
	{
		_max_sift = max_sift;
		AllocateSiftMatch();
		_have_loc[0] = _have_loc[1] = 0;
		_id_sift[0] = _id_sift[1] = -1;
		_num_sift[0] = _num_sift[1] = 1;
	}else
	{
		_max_sift = max_sift;
	}

}

void SiftMatchGL::AllocateSiftMatch()
{
	//parameters, number of sift is limited by the texture size
	if(_max_sift > GlobalUtil::_texMaxDimGL) _max_sift = GlobalUtil::_texMaxDimGL;	
	///
	int h = _max_sift / _sift_per_row; 
	int n = (GlobalUtil::_texMaxDimGL + h - 1) / GlobalUtil::_texMaxDimGL; 
	if ( n > 1) {_sift_num_stripe *= n; _sift_per_row *= n; }

	//initialize

	_texDes[0].InitTexture(_sift_per_row * _pixel_per_sift, _max_sift / _sift_per_row, 0,GL_RGBA8);
	_texDes[1].InitTexture(_sift_per_row * _pixel_per_sift, _max_sift / _sift_per_row, 0, GL_RGBA8);
	_texLoc[0].InitTexture(_sift_per_row , _max_sift / _sift_per_row, 0);
	_texLoc[1].InitTexture(_sift_per_row , _max_sift / _sift_per_row, 0);

	if(GlobalUtil::_SupportNVFloat || GlobalUtil::_SupportTextureRG)
	{
		//use single-component texture to save memory
#ifndef GL_R32F
#define GL_R32F 0x822E
#endif
		GLuint format = GlobalUtil::_SupportNVFloat ? GL_FLOAT_R_NV : GL_R32F;
		_texDot.InitTexture(_max_sift, _max_sift, 0, format);
		_texMatch[0].InitTexture(16, _max_sift / 16, 0, format);
		_texMatch[1].InitTexture(16, _max_sift / 16, 0, format);
	}else
	{
		_texDot.InitTexture(_max_sift, _max_sift, 0);
		_texMatch[0].InitTexture(16, _max_sift / 16, 0);
		_texMatch[1].InitTexture(16, _max_sift / 16, 0);
	}

}
void SiftMatchGL::InitSiftMatch()
{
	if(_initialized) return;
	GlobalUtil::InitGLParam();
	if(GlobalUtil::_GoodOpenGL == 0) return;
	AllocateSiftMatch();
#if !defined(SIFTGPU_NO_CG)
	if(_use_glsl)
	{
		LoadSiftMatchShadersGLSL();
	}
	else 
	{
		ProgramCG::InitContext();
		LoadSiftMatchShadersCG();
	}
#else
	LoadSiftMatchShadersGLSL();
#endif
	_initialized = 1; 
}


void SiftMatchGL::SetDescriptors(int index, int num, const unsigned char* descriptors, int id)
{	
	if(_initialized == 0) return;
	if (index > 1) index = 1;
	if (index < 0) index = 0;
	_have_loc[index] = 0;

	//the same feature is already set
	if(id !=-1 && id == _id_sift[index]) return ;
	_id_sift[index] = id;

	if(num > _max_sift) num = _max_sift;

	sift_buffer.resize(num * 128 /4);
	memcpy(&sift_buffer[0], descriptors, 128 * num);
	_num_sift[index] = num; 
	int w = _sift_per_row * _pixel_per_sift;
	int h = (num + _sift_per_row  - 1)/ _sift_per_row; 
	sift_buffer.resize(w * h * 4, 0);
	_texDes[index].SetImageSize(w , h);
	_texDes[index].BindTex(); 
	if(_sift_num_stripe == 1)
	{
		glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, w, h, GL_RGBA,  GL_UNSIGNED_BYTE, &sift_buffer[0]);
	}else
	{
		for(int i = 0; i < _sift_num_stripe; ++i)
		{
			int ws = _sift_per_stripe * _pixel_per_sift;
			int x = i * ws;
			int pos = i * ws * h * 4; 
			glTexSubImage2D(GlobalUtil::_texTarget, 0, x, 0, ws, h, GL_RGBA, GL_UNSIGNED_BYTE, &sift_buffer[pos]);
		}
	}
	_texDes[index].UnbindTex();

}

void SiftMatchGL::SetFeautreLocation(int index, const float* locations, int gap)
{
	if(_num_sift[index] <=0) return;
	int w = _sift_per_row ;
	int h = (_num_sift[index] + _sift_per_row  - 1)/ _sift_per_row; 
	sift_buffer.resize(_num_sift[index] * 2);
	if(gap == 0)
	{
		memcpy(&sift_buffer[0], locations, _num_sift[index] * 2 * sizeof(float));
	}else
	{
		for(int i = 0; i < _num_sift[index]; ++i)
		{
			sift_buffer[i*2] = *locations++;
			sift_buffer[i*2+1]= *locations ++;
			locations += gap;
		}
	}
	sift_buffer.resize(w * h * 2, 0);
	_texLoc[index].SetImageSize(w , h);
	_texLoc[index].BindTex(); 
	if(_sift_num_stripe == 1)
	{
		glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, w, h, GL_LUMINANCE_ALPHA , GL_FLOAT , &sift_buffer[0]);
	}else
	{
		for(int i = 0; i < _sift_num_stripe; ++i)
		{
			int ws = _sift_per_stripe;
			int x = i * ws;
			int pos = i * ws * h * 2; 
			glTexSubImage2D(GlobalUtil::_texTarget, 0, x, 0, ws, h, GL_LUMINANCE_ALPHA , GL_FLOAT, &sift_buffer[pos]);
		}
	}
	_texLoc[index].UnbindTex();
	_have_loc[index] = 1;
}

void SiftMatchGL::SetDescriptors(int index, int num, const float* descriptors, int id)
{	
	if(_initialized == 0) return;
	if (index > 1) index = 1;
	if (index < 0) index = 0;
	_have_loc[index] = 0;

	//the same feature is already set
	if(id !=-1 && id == _id_sift[index]) return ;
	_id_sift[index] = id; 

	if(num > _max_sift) num = _max_sift;

	sift_buffer.resize(num * 128 /4);
	unsigned char * pub = (unsigned char*) &sift_buffer[0];
	for(int i = 0; i < 128 * num; ++i)
	{
		pub[i] = int(512 * descriptors[i] + 0.5);
	}
	_num_sift[index] = num; 
	int w = _sift_per_row * _pixel_per_sift;
	int h = (num + _sift_per_row  - 1)/ _sift_per_row; 
	sift_buffer.resize(w * h * 4, 0);
	_texDes[index].SetImageSize(w, h);
	_texDes[index].BindTex();
	if(_sift_num_stripe == 1)
	{
		glTexSubImage2D(GlobalUtil::_texTarget, 0, 0, 0, w, h, GL_RGBA,  GL_UNSIGNED_BYTE, &sift_buffer[0]);
	}else
	{
		for(int i = 0; i < _sift_num_stripe; ++i)
		{
			int ws = _sift_per_stripe * _pixel_per_sift;
			int x = i * ws;
			int pos = i * ws * h * 4; 
			glTexSubImage2D(GlobalUtil::_texTarget, 0, x, 0, ws, h, GL_RGBA, GL_UNSIGNED_BYTE, &sift_buffer[pos]);
		}
	}
	_texDes[index].UnbindTex();
}

void SiftMatchGL::LoadSiftMatchShadersCG()
{
#if !defined(SIFTGPU_NO_CG)
	ProgramCG * program;
	char buffer[10240];
	ostrstream out(buffer, 10240);
	out <<	"#define SIFT_PER_STRIPE " << _sift_per_stripe << "\n" 
			"#define PIXEL_PER_SIFT " << _pixel_per_sift << "\n"
			"void main(uniform samplerRECT tex1,		\n"
			"uniform samplerRECT tex2,			\n"
			"uniform float2	size,  \n"
			"in		float2	pos : WPOS,	\n"
			"out	float4  result:COLOR0)		\n"
		    "{\n"
		<<	"   float4 val = float4(0, 0, 0, 0), data1, buf;\n"
			"   float2 index = pos.yx; \n"
			"   float2 stripe_size = size.xy * SIFT_PER_STRIPE;\n"
			"   float2 temp_div1 = index / stripe_size;\n"
			"   float2 stripe_index = floor(temp_div1);\n"
			"   index = floor(stripe_size * (temp_div1 - stripe_index));\n"
			"   float2 temp_div2 = index / SIFT_PER_STRIPE;\n"
			"   float2 temp_floor2 = floor(temp_div2);\n"
			"   float2 index_v = temp_floor2 + 0.5;\n "
			"   float2 index_h = SIFT_PER_STRIPE* (temp_div2 - temp_floor2);\n"
			"   float2 tx = (index_h + stripe_index * SIFT_PER_STRIPE)* PIXEL_PER_SIFT + 0.5;\n"
			"   float2 tpos1, tpos2; \n"
			"   float4 tpos = float4(tx, index_v);\n"
			//////////////////////////////////////////////////////
			"   for(int i = 0; i < PIXEL_PER_SIFT; ++i){\n"
			"       buf = texRECT(tex2, tpos.yw);\n"
			"       data1 = texRECT(tex1, tpos.xz);\n"
			"       val += data1 * buf;\n"
			"       tpos.xy = tpos.xy + float2(1.0, 1.0);\n"
			"   }\n"
			"   const float factor = 0.248050689697265625; \n"
			"   result =float4(dot(val, factor.xxxx), index,  0);\n"
			"}"
		<<	'\0';

	s_multiply = program= new ProgramCG(buffer); 

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

	out.seekp(ios::beg);
	out <<	"#define SIFT_PER_STRIPE " << _sift_per_stripe << "\n" 
			"#define PIXEL_PER_SIFT " << _pixel_per_sift << "\n"
			"void main(uniform samplerRECT tex1,		\n"
			"uniform samplerRECT tex2,			\n"
			"uniform samplerRECT texL1, \n"
			"uniform samplerRECT texL2, \n"
			"uniform float3x3 H, \n"
			"uniform float3x3 F, \n"
			"uniform float4	size,  \n"
			"in		float2	pos : WPOS,	\n"
			"out	float4  result:COLOR0)		\n"
		    "{\n"
		<<	"   float4 val = float4(0, 0, 0, 0), data1, buf;\n"
			"   float2 index = pos.yx; \n"
			//"   float2 stripe_size = size.xy * SIFT_PER_STRIPE;\n"
			//"   float2 stripe_index = floor(index /stripe_size);\n"
			//"   index = floor(fmod(index,  stripe_size));\n"
			//"   float2 index_v = floor(index / SIFT_PER_STRIPE) + 0.5;\n "
			//"   float2 index_h = fmod(index, SIFT_PER_STRIPE);\n"
			"   float2 stripe_size = size.xy * SIFT_PER_STRIPE;\n"
			"   float2 temp_div1 = index / stripe_size;\n"
			"   float2 stripe_index = floor(temp_div1);\n"
			"   index = floor(stripe_size * (temp_div1 - stripe_index));\n"
			"   float2 temp_div2 = index / SIFT_PER_STRIPE;\n"
			"   float2 temp_floor2 = floor(temp_div2);\n"
			"   float2 index_v = temp_floor2 + 0.5;\n "
			"   float2 index_h = SIFT_PER_STRIPE* (temp_div2 - temp_floor2);\n"
			//read feature location data
			"   float4 tlpos = float4((index_h + stripe_index * SIFT_PER_STRIPE) + 0.5, index_v);\n"
			"   float3 loc1 = float3(texRECT(texL1, tlpos.xz).xw, 1);\n"
			"   float3 loc2 = float3(texRECT(texL2, tlpos.yw).xw, 1);\n"
			//check the guiding homography
			"   float3 hxloc1 = mul(H, loc1);\n"
			"   float2 diff = abs(loc2.xy- hxloc1.xy/hxloc1.z);\n"
			"   float disth = max(diff.x, diff.y);\n"
			"   if(disth > size.z ) {result = float4(0, index, 0); return;}\n"

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

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

	s_guided_mult = program= new ProgramCG(buffer);

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

	//row max
	out.seekp(ios::beg);
	out <<	"#define BLOCK_WIDTH "  << 16 << "\n"
			"void main (uniform samplerRECT tex, \n"
			"uniform float3 param, in float2 pos : WPOS, \n"
			"out float4 result: COLOR0)\n"
			"{\n"
			"   float index = pos.x + floor(pos.y) * BLOCK_WIDTH; \n"
			"   float2 bestv = -1; float imax = -1;\n"
			"   for(float i = 0; i < param.x; i ++){\n "
			"		float v = texRECT(tex, float2(i + 0.5, index)).r; \n"
			"		imax = v > bestv.r ? i : imax; \n "
			"		bestv  = v > bestv.r? float2(v, bestv.r) : max(bestv, v.xx);\n "
			"	}\n"
			"	bestv = acos(min(bestv, 1.0));\n"
			"	if(bestv.x >= param.y || bestv.x >= param.z * bestv.y) imax = -1;\n"
			"	result = float4(imax, bestv, index);\n"
			"}"
		<<  '\0';
	s_row_max = program= new ProgramCG(buffer); 
	_param_rowmax_param = cgGetNamedParameter(*program, "param");

	out.seekp(ios::beg);
	out <<	"#define BLOCK_WIDTH "  << 16 << "\n"
			"void main (uniform samplerRECT tex, \n"
			"uniform float3 param, in float2 pos : WPOS, \n"
			"out float4 result: COLOR0)\n"
			"{\n"
			"	float index = pos.x + floor(pos.y) * BLOCK_WIDTH; \n"
			"	float2 bestv = -1; float imax = -1;\n"
			"	for(float i = 0; i < param.x; i ++){\n "
			"		float v = texRECT(tex, float2(index, i + 0.5)).r; \n"
			"		imax = (v > bestv.r)? i : imax; \n "
			"		bestv  = v > bestv.r? float2(v, bestv.r) : max(bestv, v.xx);\n "
			"	}\n"
			"	bestv = acos(min(bestv, 1.0));\n"
			"	if(bestv.x >= param.y || bestv.x >= param.z * bestv.y) imax = -1;\n"
			"	result = float4(imax, bestv, index);\n"
			"}"
		<<  '\0';
	s_col_max = program =new ProgramCG(buffer); 
	_param_colmax_param = cgGetNamedParameter(*program, "param");
#endif
}

void SiftMatchGL::LoadSiftMatchShadersGLSL()
{
	ProgramGLSL * program;
	char buffer[10240];
	ostrstream out(buffer, 10240);
	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; uniform vec2	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"
			"   vec2 tx = (index_h + stripe_index * vec2(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"