stim/gl/gl_spider.h

#ifndef STIM_GL_SPIDER_H
#define STIM_GL_SPIDER_H
#include <GL/glew.h>
#include <GL/glut.h>
#include <cuda.h>
#include <cuda_gl_interop.h>
#include <cudaGL.h>
#include <math.h>
#include "../gl/gl_texture.h"
#include "../visualization/camera.h"
#include "./error.h"
#include "../math/vector.h"
#include "../math/rect.h"
#include "../cuda/cost.h"
#include "../cuda/glbind.h"
#include <vector>
#include <iostream>
#include <fstream>
/* Technically since gl_spider inherits from gl_texture, we could
	call the init with a path to an image stack, and upload
	the images while creating the spider (calling init) */
namespace stim
{
template<typename T>
class gl_spider : public virtual gl_texture<T>
{
	//doen't use gl_texture really, just needs the GLuint id.
	//doesn't even need the texture iD really.
	private:
		stim::vec<float> position;  	//vector designating the position of the spider.
		stim::vec<float> direction;	//vector designating the orientation of the spider
						//always a unit vector.
		stim::vec<float> magnitude;	//magnitude of the direction vector.
						//mag[0] = length.
						//mag[1] = width.
		std::vector<stim::vec<float> > dirVectors;
		std::vector<stim::vec<float> > posVectors;
		std::vector<stim::vec<float> > magVectors;
		double currentTransform[16];
		using gl_texture<T>::texID;
		using gl_texture<T>::S;
		using gl_texture<T>::R;
		cudaArray* c_Array;
		cudaGraphicsResource_t resource;
		GLuint fboID;
		GLuint texbufferID;
		int iter;		//temporary for testing
		int numSamples;
		/// Method for finding the best scale for the spider.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. 
		void
		findOptimalDirection()
		{
			genTemplate(dirVectors, 0);
			int best = getCost();
			stim::vec<float> next = dirVectors[best];
		//	next[0] = next[0]*		
			
		}
		/// Method for finding the best direction for the spider.
		/// Not sure if necessary since the next position for the spider
		/// will be at direction * magnitude.
		void
		findOptimalPosition()
		{
			
			genTemplate(magVectors, 1);
			int best = getCost();
		}
	
		/// Method for finding the best scale for the spider.
		/// changes the x, y, z size of the spider to minimize the cost
		/// function. */
		void
		findOptimalScale()
		{
			genTemplate(magVectors, 2);
			int best = getCost();
		}
		
		void
		Optimize()
		{
			/*find the optimum direction and scale */ 
		}
		
		
		
//--------------------------------------------------------------------------//
//---------------------TEMPLATE CREATION METHODS----------------------------//
//--------------------------------------------------------------------------//
		///@param solidAngle, the size of the arc to sample.
		///Method for populating the vector arrays with sampled vectors.
		///uses the default direction vector <0,0,1>
		void
		genDirectionVectors(float solidAngle = M_PI/2)
		{
			
			vec<float> d_s = direction.cart2sph().norm();
			vec<float> temp;
			int dim = (sqrt(numSamples)-1)/2;
			float p0  	= M_PI/3;
			float dt  	= solidAngle/(2.0 * (dim + 1));
			float dp  	= p0/(2.0*(dim + 1));
	
			for(int i = -dim; i <= dim; i++){
				for(int j = -dim; j <= dim; j++){
					//Create linear index
					
					temp[0] = d_s[0]; 			//rotate vector
					temp[1] = d_s[1]+dt*i;
					temp[2] = d_s[2]+dp*j;
					
					temp = (temp.sph2cart()).norm();	//back to cart
					dirVectors.push_back(temp);
				}
			}
		}
		///@param solidAngle, the size of the arc to sample.
		///Method for populating the buffer with the sampled texture.
		///uses the default vector <0,0,0>
		void
		genPositionVectors(float delta = 0.2)
		{
			
			vec<float> temp;
			int dim = (sqrt(numSamples)-1)/2;
			stim::rect<float> samplingPlane =
				 stim::rect<float>(magnitude[0]*delta, position, direction);
			float step = 1.0/(dim);
			//Loop over the samples, keeping the original position sample
			//in the center of the resulting texture.
			for(int i = -dim; i <= dim; i++){
				for(int j = -dim; j <= dim; j++){
					//Create linear index
					temp = samplingPlane.p(
							0.5+step*i,
								 0.5+step*j
										);
					posVectors.push_back(temp);
				}
			}
		}
		///@param solidAngle, the size of the arc to sample.
		///Method for populating the buffer with the sampled texture.
		///uses the default magnitude <1,1,0>
		void
		genMagnitudeVectors(float delta = 0.5)
		{
			
			int dim = (sqrt(numSamples)-1)/2;
			float min 	= 1.0-delta;
			float max 	= 1.0+delta;
			float step	= (max-min)/(numSamples-1);
			float factor;
			vec<float> temp;
			for(int i = 0; i < numSamples; i++){
				//Create linear index
				factor = (min+step*i)*magnitude[0];
				temp[0] = factor;
				temp[1] = factor;
				magVectors.push_back(temp);	
			}
		}
		///@param vector of stim::vec in that stores all of the samplable vectors.
		///@param type, one of three operations, 0 for Direction vectors
		///	1 for Position, 2 for Magnitude.
		///Function for filling the buffer up with the data based on the vectors
		///Each vector represents a two rectangular templates.
		///Loops through all of the vectors and transfers rect. associated with it
		///Into buffer-space.
		void
		genTemplate(std::vector<stim::vec<float> > in, int type)
		{
			float x = 0.0;
			vec<float> Y(1.0,0.0,0.0);
			vec<float> pos(0.0,0.0,0.0);
			vec<float> mag(1.0, 1.0, 1.0);
                        switch(type) {
				case 0:          //Direction
					Bind();
					positionTemplate();   ///this can be placed anywhere
							      ///as long as the GL_TEXTURE
							      ///is not cleared after this call.
					for(int i = 0; i < in.size(); i++)
					{
					 	if(cos(Y.dot(in[i]))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
						else{Y[0] = 1.0; Y[1] = 0.0;}
                                         	hor = stim::rect<float>(mag,
							 pos, in[i],
                        				((Y.cross(in[i])).cross(in[i])).norm());
              				 	ver = stim::rect<float>(mag,
							 pos, in[i],
                        				hor.n());
						UpdateBuffer(x, x+i*10.0);
					 }
					 //getSample();
					 Unbind();
               				 break;
        			case 1:		 //Position
					Bind();	
					positionTemplate();
					if(cos(Y.dot(direction))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
					else{Y[0] = 1.0; Y[1] = 0.0;}
					for(int i = 0; i < in.size(); i++)
					{
                				hor = stim::rect<float>(magnitude,
							 in[i], direction,
                      					((Y.cross(direction)).cross(direction))
							.norm());
                				ver = stim::rect<float>(magnitude,
							 in[i], direction,
                        				hor.n());
						UpdateBuffer(x, x+i*10.0);
					 }
					Unbind();
                			 break;
        			case 2:		  //Scale
					Bind();	
					positionTemplate();
					if(cos(Y.dot(direction))< 0.087){                                                                             Y[0] = 0.0; Y[1] = 1.0;}
					else{Y[0] = 1.0; Y[1] = 0.0;}
                			for(int i = 0; i < in.size(); i++)
					{
						hor = stim::rect<float>(in[i],
							 position, direction,
                      					((Y.cross(direction)).cross(direction))
							.norm());
                				ver = stim::rect<float>(in[i],
							 position, direction,
                        				hor.n());
						UpdateBuffer(x, x+i*10.0);
					}
					Unbind();
					break;
        			default:
                			 std::cout << "unknown case have been passed"
                         			<< std::endl;
					 break;
			}               
		Unbind();
		CHECK_OPENGL_ERROR
		}
		///@param v_x x-coordinate in buffer-space,
		///@param v_y y-coordinate in buffer-space.
		///Samples the texturespace and places a sample in the provided coordinates
		///of bufferspace.
		void
		UpdateBuffer(float v_x, float v_y)
		{	
			float len = 10.0;
			stim::vec<float>p1; 
        	        stim::vec<float>p2; 
	                stim::vec<float>p3; 
                	stim::vec<float>p4;	
			p1 = hor.p(1,1);
			p2 = hor.p(1,0);
			p3 = hor.p(0,0);
			p4 = hor.p(0,1);
			glBegin(GL_QUADS);
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x,v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+len, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+len, v_y+len);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x, v_y+len);
			 glEnd();
			 p1 = ver.p(1,1);
			 p2 = ver.p(1,0);
			 p3 = ver.p(0,0);
			 p4 = ver.p(0,1);
		 	 glBegin(GL_QUADS);
				glTexCoord3f(
					p1[0],
					p1[1],
					p1[2]
					);
				glVertex2f(v_x+len, v_y);
				glTexCoord3f(
					p2[0],
					p2[1],
					p2[2]
					);
				glVertex2f(v_x+2*len, v_y);
				glTexCoord3f(
					p3[0],
					p3[1],
					p3[2]
					);
				glVertex2f(v_x+2*len, v_y+len);
				glTexCoord3f(
					p4[0],
					p4[1],
					p4[2]
					);
				glVertex2f(v_x+len, v_y+len);
			glEnd(); 
			CHECK_OPENGL_ERROR
		}
		
//--------------------------------------------------------------------------//
//--------------------------------GL METHODS--------------------------------//
//--------------------------------------------------------------------------//
		///@param width sets the width of the buffer.
		///@param height sets the height of the buffer.
		///Function for setting up the 2D buffer that stores the samples.
		void
		GenerateFBO(unsigned int width, unsigned int height)
		{
			glGenFramebuffers(1, &fboID);
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
			int numChannels = 1;
			unsigned char* texels = new unsigned char[width * height * numChannels];
			glGenTextures(1, &texbufferID);
			glBindTexture(GL_TEXTURE_2D, texbufferID);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_REPEAT);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
			glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
			glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE,
				 width, height, 0, GL_LUMINANCE, GL_UNSIGNED_BYTE, texels);   
			delete[] texels;
			glBindFramebuffer(GL_FRAMEBUFFER, 0); 
			glBindTexture(GL_TEXTURE_2D, 0);
		}
		///Method for controling the buffer and texture binding in order to properly
		///do the render to texture.
		void
		Bind()
		{
			float len = 10.0;
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);//set up GL buffer		
			glFramebufferTexture2D(
				GL_FRAMEBUFFER,
				GL_COLOR_ATTACHMENT0,
				GL_TEXTURE_2D,
				texbufferID,
				0);
			glBindFramebuffer(GL_FRAMEBUFFER, fboID);
			GLenum DrawBuffers[1] = {GL_COLOR_ATTACHMENT0};
			glDrawBuffers(1, DrawBuffers);
			glBindTexture(GL_TEXTURE_2D, texbufferID);
			glClearColor(1,1,1,1);
			glClear(GL_COLOR_BUFFER_BIT);
			glMatrixMode(GL_PROJECTION);
			glLoadIdentity();
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glViewport(0,0,2.0*len, numSamples*len);
			gluOrtho2D(0.0,2.0*len,0.0,numSamples*len);
			glEnable(GL_TEXTURE_3D);
			glBindTexture(GL_TEXTURE_3D, texID);
			CHECK_OPENGL_ERROR
		}
		
		///Method for Unbinding all of the texture resources
		void
		Unbind()
		{
			//Finalize GL_buffer
			glBindTexture(GL_TEXTURE_3D, 0);                      
			glDisable(GL_TEXTURE_3D);
			glBindFramebuffer(GL_FRAMEBUFFER,0);
			glBindTexture(GL_TEXTURE_2D, 0);
		}
		///Method for using the gl manipulation to alighn templates from
		///Template space (-0.5 0.5) to Texture space (0.0, 1.0),
		///Based on the position of the spider in real space (arbitrary).
		void positionTemplate()
		{
			stim::vec<float, 4> rot = getRotation(direction);
			glMatrixMode(GL_TEXTURE);
			glLoadIdentity();
			glScalef(1.0/512.0/0.6, 1.0/512.0/0.6, 1.0/426.0/2.0);
		glGetDoublev(GL_TEXTURE_MATRIX, currentTransform);
		printTransform();
			glTranslatef(position[0],
				     position[1],
				     position[2]);
		glGetDoublev(GL_TEXTURE_MATRIX, currentTransform);
		printTransform();
			glScalef(magnitude[0],
				 magnitude[1],
				 magnitude[0]);
		glGetDoublev(GL_TEXTURE_MATRIX, currentTransform);
		printTransform();
			glRotatef(rot[0], rot[1], rot[2], rot[3]);
		glGetDoublev(GL_TEXTURE_MATRIX, currentTransform);
		printTransform();
			glGetDoublev(GL_TEXTURE_MATRIX, currentTransform);
			CHECK_OPENGL_ERROR
			glMatrixMode(GL_MODELVIEW);
		}
		
		
//--------------------------------------------------------------------------//
//--------------------------------CUDA METHODS------------------------------//
//--------------------------------------------------------------------------//
		
		/// Method for registering the texture with Cuda for shared
		///	access.
		void
		createResource()
		{
			HANDLE_ERROR(
				cudaGraphicsGLRegisterImage(
					 &resource,
				 	texbufferID,
				 	GL_TEXTURE_2D,
				 	//CU_GRAPHICS_REGISTER_FLAGS_NONE)
					cudaGraphicsMapFlagsReadOnly)
			);
		} 
		
		///Method for freeing the texture from Cuda for gl access.
		void
		destroyResource()
		{
			HANDLE_ERROR(
				cudaGraphicsUnregisterResource(resource)
			);		
		}
		///Entry-point into the cuda code for calculating the cost
		///	of a given samples array (in texture form) 
		int
		getCost()
		{
			createResource();
			int cost = 	get_cost(resource, iter, numSamples);
			destroyResource();
			iter++;
			return cost;
		}
	public:
		stim::rect<float> hor;
		stim::rect<float> ver;	
//--------------------------------------------------------------------------//
//-----------------------------CONSTRUCTORS---------------------------------//
//--------------------------------------------------------------------------//
		///Default Constructor
		gl_spider
		()
		{
			setPosition(0.0,0.0,0.0);
			setDirection(0.0,0.0,1.0);
			setMagnitude(1.0);
			//numSamples = 1089;
			numSamples = 9;
		}
		///@param samples, the number of samples this spider is going to use.
		///best results if samples is can create a perfect root.
		///Default Constructor
		gl_spider
		(int samples)
		{
			setPosition(0.0,0.0,0.0);
			setDirection(0.0,0.0,1.0);
			setMagnitude(1.0);
			numSamples = samples;
		}
		///temporary constructor for convenience, will be removed in further updates.	
		gl_spider
		(float pos_x, float pos_y, float pos_z, float dir_x, float dir_y, float dir_z,
			float mag_x)
		{
			setPosition(pos_x, pos_y, pos_z);
			setDirection(dir_x, dir_y, dir_z);
			setMagnitude(mag_x);
		
		}
		///@param GLuint id texture that is going to be sampled.
		///Attached the spider to the texture with the given GLuint ID.
		///Samples in the default direction acting as the init method.
		///Also acts an init.	
		void
		attachSpider(GLuint id)
		{
			texID = id;
			iter = 0; ///for debugging purposes
			GenerateFBO(20, numSamples*10);
			genDirectionVectors();
			genPositionVectors();
			genMagnitudeVectors();
			gl_texture<T>::setDims(0.6, 0.6, 2.0);
			genTemplate(dirVectors, 0);
		}
		
//--------------------------------------------------------------------------//
//-----------------------------ACCESS METHODS-------------------------------//
//--------------------------------------------------------------------------//
		///Returns the position vector.
		vec<float>
		getPosition()
		{
			return position;
		}
	
		///Returns the direction vector.
		vec<float>
		getDirection()
		{
			return direction;
		}
		///Returns the magnitude vector.
		vec<float>
		getMagnitude()
		{
			return magnitude;
		}
	
		///@param vector pos, the new position.
		///Sets the position vector to input vector pos.
		void
		setPosition(vec<float> pos)
		{
			position = pos;
		}
		///@param x x-coordinate.
		///@param y y-coordinate.
		///@param z z-coordinate.
		///Sets the position vector to the input float coordinates x,y,z.
		void
		setPosition(float x, float y, float z)
		{
			position[0] = x;
			position[1] = y;
			position[2] = z;
		}
		
		///@param vector dir, the new direction.
		///Sets the direction vector to input vector dir.
		void
		setDirection(vec<float> dir)
		{
			direction = dir;
		}
		
		///@param x x-coordinate.
		///@param y y-coordinate.
		///@param z z-coordinate.
		///Sets the direction vector to the input float coordinates x,y,z.
		void
		setDirection(float x, float y, float z)
		{
			direction[0] = x;
			direction[1] = y;
			direction[2] = z;
		}
			
		///@param vector dir, the new direction.
		///Sets the magnitude vector to the input vector mag.	
		void
		setMagnitude(vec<float> mag)
		{
			magnitude[0] = mag[0];
			magnitude[1] = mag[0];
		}
		
		///@param mag size of the sampled region.
		///Sets the magnitude vector to the input mag for both templates.
		void
		setMagnitude(float mag)
		{
			magnitude[0] = mag;
			magnitude[1] = mag;
		}
		
		void
		setSize(int x, int y, int z)
		{
			R[1] = x;
			R[2] = y;
			R[3] = z;
		}
		
		///@param dir, the vector to which we are rotating
		///given a vector to align to, finds the required
		///axis and angle for glRotatef
		stim::vec<float, 4>
		getRotation(stim::vec<float> dir)
		{
			stim::vec<float, 4> out;
			stim::vec<float> from(0,0,1);
			out[0] = acos(from.dot(dir))*M_PI/180;
			if(out[0] < 0.0001){
				out[0] = 0.0;
				out[1] = 0.0;
				out[2] = 0.0;
				out[3] = 1.0;
			} else {
				stim::vec<float> temp;
				temp = (from.cross(dir)).norm();
				out[1] = temp[0];
				out[2] = temp[1];
				out[3] = temp[2];
			}
			return out;
		}
		
		///temporary method for visualization.
		GLuint
		getFB()
		{
			return fboID;
		}
//--------------------------------------------------------------------------//
//-----------------------------TEMPORARY METHODS----------------------------//
//--------------------------------------------------------------------------//
		///temporary Method necessary for visualization and testing.
		void
		Update()
		{
			vec<float> Y(1.0,0.0,0.0);
			if(cos(Y.dot(direction))< 0.087){
				Y[0] = 0.0; Y[1] = 1.0;}
			hor = stim::rect<float>(magnitude, position, direction.norm(),
				((Y.cross(direction)).cross(direction)).norm());
			ver = stim::rect<float>(magnitude, position, direction.norm(),
				 hor.n());
		}
		void
		Step()
		{
			findOptimalDirection();
		}
		void
		printTransform()
		{
		for(int i = 0; i < 4; i++){
			std::cout << "[" << currentTransform[i] << "]" <<
				     "[" << currentTransform[i+4] << "]" <<
				     "[" << currentTransform[i+8] << "]" <<
				     "[" << currentTransform[i+12] << "]" <<
			std::endl;
				}
		}
		/* Method for initializing the cuda devices, necessary only
			there are multiple cuda devices */
		void
		initCuda()
		{	
			stim::cudaSetDevice();
			//GLint max;
			//glGetIntegerv(GL_MAX_TEXTURE_SIZE, &max);
			//std::cout << max << std::endl;
		}
};
}
#endif