programglsl.cpp

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

		out<<"coord = gl_TexCoord[0].xy + vec2(float("<<i-nhpixel<<"),0);\n";
		out<<"pc=texture2DRect(tex, coord);\n";
		if(GlobalUtil::_PreciseBorder)		out<<"if(coord.x < 0) pc = pc.rrbb;\n";
		//for each sub-pixel j  in center, the weight of sub-pixel k 
		xw = (i - nhpixel)*2;
		for( j = 0; j < 3; j++)
		{
			xwn = xw  + j  -1;
			weight[j] = xwn < -halfwidth || xwn > halfwidth? 0 : pf[xwn];
		}
		if(weight[1] == 0.0)
		{
			out<<"result += vec4("<<weight[2]<<","<<weight[0]<<","<<weight[2]<<","<<weight[0]<<")*pc.grab;\n";
		}
		else
		{
			out<<"result += vec4("<<weight[1]<<", "<<weight[0]<<", "<<weight[1]<<", "<<weight[0]<<")*pc.rrbb;\n";
			out<<"result += vec4("<<weight[2]<<", "<<weight[1]<<", "<<weight[2]<<", "<<weight[1]<<")*pc.ggaa;\n";
		}	
	
	}
	out<<"gl_FragColor = result;}\n"<<'\0';
	
	return new ProgramGLSL( buffer);


}


ProgramGPU* FilterGLSL::CreateFilterVPK(float kernel[], float offset[], int height)
{

	//both h and v are packed...
	int i, j, yw, ywn;

	int halfh  = height >>1;
	float * pf = kernel + halfh;
	int nhpixel = (halfh+1)>>1;	//how many neighbour pixels need to be looked up
	int npixel  = (nhpixel<<1)+1;//
	char buffer[10240];
	float weight[3];
	ostrstream out(buffer, 10240);
	out<<setprecision(8);

	out<<  "uniform sampler2DRect tex;";
	out<< "\nvoid main(void){ vec4 result = vec4(0, 0, 0, 0);\n";
	///use multi texture coordinate because nhpixels can be at most 3
	out<<"vec4 pc; vec2 coord;\n";
	for( i = 0 ; i < npixel ; i++)
	{
		out<<"coord = gl_TexCoord[0].xy + vec2(0, float("<<i-nhpixel<<"));\n";
		out<<"pc=texture2DRect(tex, coord);\n";
		if(GlobalUtil::_PreciseBorder)	out<<"if(coord.y < 0) pc = pc.rgrg;\n";

		//for each sub-pixel j  in center, the weight of sub-pixel k 
		yw = (i - nhpixel)*2;
		for( j = 0; j < 3; j++)
		{
			ywn = yw + j  -1;
			weight[j] = ywn < -halfh || ywn > halfh? 0 : pf[ywn];
		}
		if(weight[1] == 0.0)
		{
			out<<"result += vec4("<<weight[2]<<","<<weight[2]<<","<<weight[0]<<","<<weight[0]<<")*pc.barg;\n";
		}else
		{
			out<<"result += vec4("<<weight[1]<<","<<weight[1]<<","<<weight[0]<<","<<weight[0]<<")*pc.rgrg;\n";
			out<<"result += vec4("<<weight[2]<<","<<weight[2]<<","<<weight[1]<<","<<weight[1]<<")*pc.baba;\n";
		}

	}
	out<<"gl_FragColor = result;}\n"<<'\0';

	return new ProgramGLSL( buffer);
}








void ShaderBagGLSL::LoadFixedShaders()
{


	s_gray = new ProgramGLSL( 
		"uniform sampler2DRect rgbTex; void main(void){\n"
		"float intensity = dot(vec3(0.299, 0.587, 0.114), texture2DRect(rgbTex,gl_TexCoord[0].st ).rgb);\n"
		"gl_FragColor = vec4(intensity, intensity, intensity, 1.0);}");


	s_debug = new ProgramGLSL( "void main(void){gl_FragColor.rg =  gl_TexCoord[0].st;}");


	s_sampling = new ProgramGLSL(
		"uniform sampler2DRect tex; void main(void){gl_FragColor.rg= texture2DRect(tex, gl_TexCoord[0].st).rg;}");

	//
	s_grad_pass = new ProgramGLSL(
	"uniform sampler2DRect tex; void main ()\n"
	"{\n"
	"	vec4 v1, v2, gg;\n"
	"	vec4 cc  = texture2DRect(tex, gl_TexCoord[0].xy);\n"
	"	gg.x = texture2DRect(tex, gl_TexCoord[1].xy).r;\n"
	"	gg.y = texture2DRect(tex, gl_TexCoord[2].xy).r;\n"
	"	gg.z = texture2DRect(tex, gl_TexCoord[3].xy).r;\n"
	"	gg.w = texture2DRect(tex, gl_TexCoord[4].xy).r;\n"
	"	vec2 dxdy = (gg.yw - gg.xz); \n"
	"	float grad = 0.5*length(dxdy);\n"
	"	float theta = grad==0? 0.0: atan(dxdy.y, dxdy.x);\n"
	"	gl_FragData[0] = vec4(cc.rg, grad, theta);\n"
	"}\n\0");

	ProgramGLSL * program;
	s_margin_copy = program = new ProgramGLSL(
	"uniform sampler2DRect tex; uniform vec2 truncate;\n"
	"void main(){ gl_FragColor = texture2DRect(tex, min(gl_TexCoord[0].xy, truncate)); }");

	_param_margin_copy_truncate = glGetUniformLocation(*program, "truncate");


	GlobalUtil::_OrientationPack2 = 0;
	LoadOrientationShader();

	if(s_orientation == NULL)
	{
		//Load a simplified version if the right version is not supported
		s_orientation = program =  new ProgramGLSL(
		"uniform sampler2DRect fTex; uniform sampler2DRect oTex;\n"
	"	uniform float size; void main(){\n"
	"	vec4 cc = texture2DRect(fTex, gl_TexCoord[0].st);\n"
	"	vec4 oo = texture2DRect(oTex, cc.rg);\n"
	"	gl_FragColor.rg = cc.rg;\n"
	"	gl_FragColor.b = oo.a;\n"
	"	gl_FragColor.a = size;}");  

		_param_orientation_gtex = glGetUniformLocation(*program, "oTex");
		_param_orientation_size = glGetUniformLocation(*program, "size");

		GlobalUtil::_MaxOrientation = 0;
		GlobalUtil::_FullSupported = 0;
		std::cerr<<"Orientation simplified on this hardware"<<endl;
	}


	if(GlobalUtil::_DescriptorPPT) LoadDescriptorShader();
	if(s_descriptor_fp == NULL) 
	{
		GlobalUtil::_DescriptorPPT = GlobalUtil::_FullSupported = 0; 
		std::cerr<<"Descriptor ignored on this hardware"<<endl;
	}

	s_zero_pass = new ProgramGLSL("void main(){gl_FragColor = vec4(0.0);}");
}


void ShaderBagGLSL::LoadDisplayShaders()
{
	s_copy_key = new ProgramGLSL(
		"uniform sampler2DRect tex; void main(){\n"
	"gl_FragColor.rg= texture2DRect(tex, gl_TexCoord[0].st).rg; gl_FragColor.ba = vec2(0.0,1.0);	}");
	

	ProgramGLSL * program;
	s_vertex_list = program = new ProgramGLSL(
	"uniform vec4 sizes; uniform sampler2DRect tex;\n"
	"void main(void){\n"
	"float fwidth = sizes.y; float twidth = sizes.z; float rwidth = sizes.w; \n"
	"float index = 0.1*(fwidth*floor(gl_TexCoord[0].y) + gl_TexCoord[0].x);\n"
	"float px = mod(index, twidth);\n"
	"vec2 tpos= floor(vec2(px, index*rwidth))+0.5;\n"
	"vec4 cc = texture2DRect(tex, tpos );\n"
	"float size = 3.0f * cc.a; //sizes.x;// \n"
	"gl_FragColor.zw = vec2(0.0, 1.0);\n"
	"if(cc.x<=0 || cc.y <=0) {gl_FragColor.xy = cc.xy; }\n"
	"else {float type = fract(px);\n"
	"vec2 dxy; \n"
	"dxy.x = type < 0.1 ? 0.0 : ((type <0.5 || type > 0.9)? size : -size);\n"
	"dxy.y = type < 0.2 ? 0.0 : ((type < 0.3 || type > 0.7 )? -size :size); \n"
	"float s = sin(cc.b); float c = cos(cc.b); \n"
	"gl_FragColor.x = cc.x + c*dxy.x-s*dxy.y;\n"
	"gl_FragColor.y = cc.y + c*dxy.y+s*dxy.x;}\n}\n");


	_param_genvbo_size = glGetUniformLocation(*program, "sizes");
	s_display_gaussian =  new ProgramGLSL(
	"uniform sampler2DRect tex; void main(void){float r = texture2DRect(tex, gl_TexCoord[0].st).r;\n"
	"gl_FragColor = vec4(r, r, r, 1);}" );

	s_display_dog =  new ProgramGLSL(
	"uniform sampler2DRect tex; void main(void){float g = 0.5+(20.0*texture2DRect(tex, gl_TexCoord[0].st).g);\n"
	"gl_FragColor = vec4(g, g, g, 0.0);}" );

	s_display_grad = new ProgramGLSL(
		"uniform sampler2DRect tex; void main(void){\n"
    "	vec4 cc = texture2DRect(tex, gl_TexCoord[0].st);gl_FragColor = vec4(5.0* cc.bbb, 1.0);}");

	s_display_keys= new ProgramGLSL(
		"uniform sampler2DRect tex; void main(void){\n"
	"	vec4 cc = texture2DRect(tex, gl_TexCoord[0].st);\n"
	"	if(cc.r ==0.0) discard; gl_FragColor =  (cc.r==1.0? vec4(1.0, 0.0, 0,1.0):vec4(0.0,1.0,0.0,1.0));}");
}

void ShaderBagGLSL::LoadKeypointShader(float threshold, float edge_threshold)
{

	char buffer[10240];
	float threshold0 = threshold* (GlobalUtil::_SubpixelLocalization?0.8f:1.0f);	
	float threshold1 = threshold;
	float threshold2 = (edge_threshold+1)*(edge_threshold+1)/edge_threshold;
	ostrstream out(buffer, 10240);
	streampos pos;

	//tex(X)(Y)
	//X: (CLR) (CENTER 0, LEFT -1, RIGHT +1)  
	//Y: (CDU) (CENTER 0, DOWN -1, UP    +1) 

	out <<	"#define THRESHOLD0 " << threshold0 << "\n"
			"#define THRESHOLD1 " << threshold1 << "\n"
			"#define THRESHOLD2 " << threshold2 << "\n";

	out<<
	"uniform sampler2DRect tex, texU, texD; void main ()\n"
	"{\n"
	"	vec4 v1, v2, gg, temp;\n"
	"	vec2 TexRU = vec2(gl_TexCoord[2].x, gl_TexCoord[4].y); \n"
	"	vec4 cc  = texture2DRect(tex, gl_TexCoord[0].xy);\n"
	"	temp =  texture2DRect(tex, gl_TexCoord[1].xy);\n"
	"	v1.x =  temp.g;			gg.x = temp.r;\n"
	"	temp = texture2DRect(tex, gl_TexCoord[2].xy) ;\n"
	"	v1.y = temp.g;			gg.y = temp.r;\n"
	"	temp = texture2DRect(tex, gl_TexCoord[3].xy) ;\n"
	"	v1.z = temp.g;			gg.z = temp.r;\n"
	"	temp = texture2DRect(tex, gl_TexCoord[4].xy) ;\n"
	"	v1.w = temp.g;			gg.w = temp.r;\n"
	"	v2.x = texture2DRect(tex, gl_TexCoord[5].xy).g;\n"
	"	v2.y = texture2DRect(tex, gl_TexCoord[6].xy).g;\n"
	"	v2.z = texture2DRect(tex, gl_TexCoord[7].xy).g;\n"
	"	v2.w = texture2DRect(tex, TexRU.xy).g;\n"
	"	vec2 dxdy = (gg.yw - gg.xz); \n"
	"	float grad = 0.5*length(dxdy);\n"
	"	float theta = grad==0? 0.0: atan(dxdy.y, dxdy.x);\n"
	"	gl_FragData[0] = vec4(cc.rg, grad, theta);\n"

	//test against 8 neighbours
	//use variable to identify type of extremum
	//1.0 for local maximum and 0.5 for minimum
	<<
	"	float dog = 0.0; \n"
	"	gl_FragData[1] = vec4(0, 0, 0, 0); \n"
	"	dog = cc.g > THRESHOLD0 && all(greaterThan(cc.gggg, max(v1, v2)))?1.0: 0.0;\n"
	"	dog = cc.g < -THRESHOLD0 && all(lessThan(cc.gggg, min(v1, v2)))?0.5: dog;\n"
	"	if(dog == 0.0) return;\n";

	pos = out.tellp();
	//do edge supression first.. 
	//vector v1 is < (-1, 0), (1, 0), (0,-1), (0, 1)>
	//vector v2 is < (-1,-1), (-1,1), (1,-1), (1, 1)>

	out<<
	"	float fxx, fyy, fxy; \n"
	"	vec4 D2 = v1.xyzw - cc.gggg;\n"
	"	vec2 D4 = v2.xw - v2.yz;\n"
	"	fxx = D2.x + D2.y;\n"
	"	fyy = D2.z + D2.w;\n"
	"	fxy = 0.25*(D4.x + D4.y);\n"
	"	float fxx_plus_fyy = fxx + fyy;\n"
	"	float score_up = fxx_plus_fyy*fxx_plus_fyy; \n"
	"	float score_down = (fxx*fyy - fxy*fxy);\n"
	"	if( score_down <= 0 || score_up > THRESHOLD2 * score_down)return;\n";

	//...
	out<<" \n"
	"	vec2 D5 = 0.5*(v1.yw-v1.xz); \n"
	"	float fx = D5.x, fy = D5.y ; \n"
	"	float fs, fss , fxs, fys ; \n"
	"	vec2 v3; vec4 v4, v5, v6;\n"
	//read 9 pixels of upper level
	<<
	"	v3.x = texture2DRect(texU, gl_TexCoord[0].xy).g;\n"
	"	v4.x = texture2DRect(texU, gl_TexCoord[1].xy).g;\n"
	"	v4.y = texture2DRect(texU, gl_TexCoord[2].xy).g;\n"
	"	v4.z = texture2DRect(texU, gl_TexCoord[3].xy).g;\n"
	"	v4.w = texture2DRect(texU, gl_TexCoord[4].xy).g;\n"
	"	v6.x = texture2DRect(texU, gl_TexCoord[5].xy).g;\n"
	"	v6.y = texture2DRect(texU, gl_TexCoord[6].xy).g;\n"
	"	v6.z = texture2DRect(texU, gl_TexCoord[7].xy).g;\n"
	"	v6.w = texture2DRect(texU, TexRU.xy).g;\n"
	//compare with 9 pixels of upper level
	//read and compare with 9 pixels of lower level
	//the maximum case
	<<
	"	if(dog == 1.0)\n"
	"	{\n"
	"		if(cc.g < v3.x || any(lessThan(cc.gggg, v4)) ||any(lessThan(cc.gggg, v6)))return; \n"
	"		v3.y = texture2DRect(texD, gl_TexCoord[0].xy).g;\n"
	"		v5.x = texture2DRect(texD, gl_TexCoord[1].xy).g;\n"
	"		v5.y = texture2DRect(texD, gl_TexCoord[2].xy).g;\n"
	"		v5.z = texture2DRect(texD, gl_TexCoord[3].xy).g;\n"
	"		v5.w = texture2DRect(texD, gl_TexCoord[4].xy).g;\n"
	"		v6.x = texture2DRect(texD, gl_TexCoord[5].xy).g;\n"
	"		v6.y = texture2DRect(texD, gl_TexCoord[6].xy).g;\n"
	"		v6.z = texture2DRect(texD, gl_TexCoord[7].xy).g;\n"
	"		v6.w = texture2DRect(texD, TexRU.xy).g;\n"
	"		if(cc.g < v3.y || any(lessThan(cc.gggg, v5)) ||any(lessThan(cc.gggg, v6)))return; \n"
	"	}\n"
	//the minimum case
	<<
	"	else{\n"
	"	if(cc.g > v3.x || any(greaterThan(cc.gggg, v4)) ||any(greaterThan(cc.gggg, v6)))return; \n"
	"		v3.y = texture2DRect(texD, gl_TexCoord[0].xy).g;\n"
	"		v5.x = texture2DRect(texD, gl_TexCoord[1].xy).g;\n"
	"		v5.y = texture2DRect(texD, gl_TexCoord[2].xy).g;\n"
	"		v5.z = texture2DRect(texD, gl_TexCoord[3].xy).g;\n"
	"		v5.w = texture2DRect(texD, gl_TexCoord[4].xy).g;\n"
	"		v6.x = texture2DRect(texD, gl_TexCoord[5].xy).g;\n"
	"		v6.y = texture2DRect(texD, gl_TexCoord[6].xy).g;\n"
	"		v6.z = texture2DRect(texD, gl_TexCoord[7].xy).g;\n"
	"		v6.w = texture2DRect(texD, TexRU.xy).g;\n"
	"		if(cc.g > v3.y || any(greaterThan(cc.gggg, v5)) ||any(greaterThan(cc.gggg, v6)))return; \n"
	"	}\n";

	if(GlobalUtil::_SubpixelLocalization)

	// sub-pixel localization FragData1 = vec4(dog, 0, 0, 0); return;
	out <<
	"	fs = 0.5*( v3.x - v3.y );  \n"
	"	fss = v3.x + v3.y - cc.g - cc.g;\n"
	"	fxs = 0.25 * ( v4.y + v5.x - v4.x - v5.y);\n"
	"	fys = 0.25 * ( v4.w + v5.z - v4.z - v5.w);\n"
	
	// 
	// let dog difference be quatratic function  of dx, dy, ds; 
	// df(dx, dy, ds) = fx * dx + fy*dy + fs * ds + 
	//				  + 0.5 * ( fxx * dx * dx + fyy * dy * dy + fss * ds * ds)
	//				  + (fxy * dx * dy + fxs * dx * ds + fys * dy * ds)
	// (fx, fy, fs, fxx, fyy, fss, fxy, fxs, fys are the derivatives)
	
	//the local extremum satisfies
	// df/dx = 0, df/dy = 0, df/dz = 0
	
	//that is 
	// |-fx|     | fxx fxy fxs |   |dx|
	// |-fy|  =  | fxy fyy fys | * |dy|
	// |-fs|     | fxs fys fss |   |ds|
	// need to solve dx, dy, ds

	// Use Gauss elimination to solve the linear system
    <<
	"	vec3 dxys = vec3(0.0);		\n"
	"	vec4 A0, A1, A2 ;			\n"
	"	A0 = vec4(fxx, fxy, fxs, -fx);	\n"
	"	A1 = vec4(fxy, fyy, fys, -fy);	\n"
	"	A2 = vec4(fxs, fys, fss, -fs);	\n"
	"	vec3 x3 = abs(vec3(fxx, fxy, fxs));		\n"
	"	float maxa = max(max(x3.x, x3.y), x3.z);	\n"
	"	if(maxa >= 1e-10 ) {						\n"
	"		if(x3.y ==maxa )							\n"
	"		{											\n"
	"			vec4 TEMP = A1; A1 = A0; A0 = TEMP;	\n"
	"		}else if( x3.z == maxa )					\n"
	"		{											\n"
	"			vec4 TEMP = A2; A2 = A0; A0 = TEMP;	\n"
	"		}											\n"
	"		A0 /= A0.x;									\n"
	"		A1 -= A1.x * A0;							\n"
	"		A2 -= A2.x * A0;							\n"
	"		vec2 x2 = abs(vec2(A1.y, A2.y));		\n"
	"		if( x2.y > x2.x )							\n"
	"		{											\n"
	"			vec3 TEMP = A2.yzw;					\n"
	"			A2.yzw = A1.yzw;						\n"
	"			A1.yzw = TEMP;							\n"
	"			x2.x = x2.y;							\n"
	"		}											\n"
	"		if(x2.x >= 1e-10) {						\n"
	"			A1.yzw /= A1.y;								\n"
	"			A2.yzw -= A2.y * A1.yzw;					\n"
	"			if(abs(A2.z) >= 1e-10) {		\n"
	// compute dx, dy, ds: 
	<<
	"				\n"
	"				dxys.z = A2.w /A2.z;				    \n"
	"				dxys.y = A1.w - dxys.z*A1.z;			    \n"
	"				dxys.x = A0.w - dxys.z*A0.z - dxys.y*A0.y;	\n"

	//one more threshold which I forgot in versions prior to 286
	<<
	"				bool dog_test = (abs(cc.g + 0.5*dot(vec3(fx, fy, fs), dxys ))<= THRESHOLD1) ;\n"
	"				if(dog_test || any(greaterThan(abs(dxys), vec3(1.0)))) dog = 0;\n"
	"			}\n"
	"		}\n"
	"	}\n"
    //keep the point when the offset is less than 1
	<<
	"	gl_FragData[1] = vec4( dog, dxys); \n";
	else

	out<<
	"	gl_FragData[1] =  vec4( dog, 0, 0, 0) ;	\n";

	out<<