From 7006df5f3b680618a1597310683f8c334ca527fb Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Mon, 12 May 2014 13:05:19 -0500 Subject: [PATCH] reformat of directory structure --- rts/biology/fibernet.h | 0 rts/cuda/glbind.h | 77 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/cuda/memory.h | 2 ++ rts/cuda/threads.h | 34 ++++++++++++++++++++++++++++++++++ rts/envi/envi.h | 2 +- rts/envi/envi_header.h | 10 +++------- rts/gl/error.h | 15 +++++++++++++++ rts/gl/rtsSourceCode.h | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/gl/rts_glShaderObject.h | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/gl/rts_glShaderProgram.h | 295 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/gl/rts_glShaderUniform.h | 105 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/gl/rts_glUtilities.h | 12 ++++++++++++ rts/gl/texture.h | 277 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/math/complex.h | 152 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------- rts/math/function.h | 123 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/math/matrix.h | 18 +++++++++++------- rts/math/point.h | 154 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------- rts/math/quad.h | 105 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------- rts/math/quaternion.h | 109 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------- rts/math/spherical_bessel.h | 10 +++++----- rts/math/triangle.h | 189 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/math/vector.h | 72 ++++++++++++++++++++++++++++++++++++++++++------------------------------ rts/optics/material.h | 107 ++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------- rts/source/colormap.cpp | 37 +++++++++++++++++++++++++++++++++++++ rts/visualization/camera.h | 198 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rts/visualization/colormap.h | 321 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 26 files changed, 2254 insertions(+), 350 deletions(-) create mode 100644 rts/biology/fibernet.h create mode 100644 rts/cuda/glbind.h create mode 100644 rts/cuda/memory.h create mode 100644 rts/cuda/threads.h create mode 100755 rts/gl/error.h create mode 100755 rts/gl/rtsSourceCode.h create mode 100755 rts/gl/rts_glShaderObject.h create mode 100755 rts/gl/rts_glShaderProgram.h create mode 100755 rts/gl/rts_glShaderUniform.h create mode 100755 rts/gl/rts_glUtilities.h create mode 100755 rts/gl/texture.h create mode 100644 rts/math/function.h create mode 100644 rts/math/triangle.h create mode 100644 rts/source/colormap.cpp create mode 100755 rts/visualization/camera.h create mode 100644 rts/visualization/colormap.h diff --git a/rts/biology/fibernet.h b/rts/biology/fibernet.h new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/rts/biology/fibernet.h diff --git a/rts/cuda/glbind.h b/rts/cuda/glbind.h new file mode 100644 index 0000000..483c504 --- /dev/null +++ b/rts/cuda/glbind.h @@ -0,0 +1,77 @@ +#ifndef RTS_GL_BIND_H +#define RTS_GL_BIND_H + +#include +#include + +#include +#include + +#include +#include "cuda_gl_interop.h" +#include "rts/gl/error.h" + + +static void rtsInitGLEW() +{ + //Initialize the GLEW toolkit + + GLenum err = glewInit(); + if(GLEW_OK != err) + { + printf("Error starting GLEW."); + } + fprintf(stdout, "Status: Using GLEW %s\n", glewGetString(GLEW_VERSION)); +} + +static void rts_cudaSetDevice(int major = 1, int minor = 3) +{ + cudaDeviceProp prop; + int dev; + + //find a CUDA device that can handle an offscreen buffer + int num_gpu; + HANDLE_ERROR(cudaGetDeviceCount(&num_gpu)); + printf("Number of CUDA devices detected: %d\n", num_gpu); + memset(&prop, 0, sizeof(cudaDeviceProp)); + prop.major=major; + prop.minor=minor; + HANDLE_ERROR(cudaChooseDevice(&dev, &prop)); + HANDLE_ERROR(cudaGetDeviceProperties(&prop, dev)); + HANDLE_ERROR(cudaGLSetGLDevice(dev)); +} + +static void* rts_cudaMapResource(cudaGraphicsResource* cudaBufferResource) +{ + //this function takes a predefined CUDA resource and maps it to a pointer + void* buffer; + HANDLE_ERROR(cudaGraphicsMapResources(1, &cudaBufferResource, NULL)); + size_t size; + HANDLE_ERROR(cudaGraphicsResourceGetMappedPointer( (void**)&buffer, &size, cudaBufferResource)); + return buffer; +} +static void rts_cudaUnmapResource(cudaGraphicsResource* resource) +{ + //this function unmaps the CUDA resource so it can be used by OpenGL + HANDLE_ERROR(cudaGraphicsUnmapResources(1, &resource, NULL)); +} + +static void rts_cudaCreateRenderBuffer(GLuint &glBufferName, cudaGraphicsResource* &cudaBufferResource, int resX, int resY) +{ + //delete the previous buffer name and resource + if(cudaBufferResource != 0) + HANDLE_ERROR(cudaGraphicsUnregisterResource(cudaBufferResource)); + if(glBufferName != 0) + glDeleteBuffers(1, &glBufferName); + + //generate an OpenGL offscreen buffer + glGenBuffers(1, &glBufferName); + + //bind the buffer - directs all calls to this buffer + glBindBuffer(GL_PIXEL_UNPACK_BUFFER, glBufferName); + glBufferData(GL_PIXEL_UNPACK_BUFFER, resX * resY * sizeof(uchar3), NULL, GL_DYNAMIC_DRAW_ARB); + CHECK_OPENGL_ERROR + HANDLE_ERROR(cudaGraphicsGLRegisterBuffer(&cudaBufferResource, glBufferName, cudaGraphicsMapFlagsNone)); +} + +#endif diff --git a/rts/cuda/memory.h b/rts/cuda/memory.h new file mode 100644 index 0000000..4cb183c --- /dev/null +++ b/rts/cuda/memory.h @@ -0,0 +1,2 @@ +#include + diff --git a/rts/cuda/threads.h b/rts/cuda/threads.h new file mode 100644 index 0000000..14fa22b --- /dev/null +++ b/rts/cuda/threads.h @@ -0,0 +1,34 @@ +#include "cuda_runtime.h" +#include "device_launch_parameters.h" +#include "rts/cuda/callable.h" + +#ifndef CUDA_THREADS_H +#define CUDA_THREADS_H + +#define MAX_GRID 65535 + +__device__ unsigned int ThreadIndex1D() +{ + return blockIdx.y * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x; +} + +dim3 GenGrid1D(unsigned int N, unsigned int blocksize = 128) +{ + dim3 dimgrid; + + dimgrid.x = (N + blocksize - 1)/blocksize; + dimgrid.y = 1; + dimgrid.z = 1; + + if(dimgrid.x > MAX_GRID) + { + dimgrid.y = (dimgrid.x + MAX_GRID - 1) / MAX_GRID; + dimgrid.x = MAX_GRID; + } + + return dimgrid; + +} + + +#endif diff --git a/rts/envi/envi.h b/rts/envi/envi.h index 1f0d09e..df82b8f 100644 --- a/rts/envi/envi.h +++ b/rts/envi/envi.h @@ -173,7 +173,7 @@ class EnviFile exit(1); } - float r, v0, v1; + float v0, v1; for(int n=0; n result; - - double fentry; - std::string entry; size_t i; do { i = sequence.find_first_of(','); entry = sequence.substr(0, i); - fentry = atof(entry.c_str()); sequence = sequence.substr(i+1); result.push_back(atof(entry.c_str())); //std::cout< 0) { outfile<<"band names = {"< +#include +#include + +#define CHECK_OPENGL_ERROR \ +{ GLenum error; \ + while ( (error = glGetError()) != GL_NO_ERROR) { \ + printf( "OpenGL ERROR: %s\nCHECK POINT: %s (line %d)\n", gluErrorString(error), __FILE__, __LINE__ ); \ + } \ +} + +#endif \ No newline at end of file diff --git a/rts/gl/rtsSourceCode.h b/rts/gl/rtsSourceCode.h new file mode 100755 index 0000000..0ff6c25 --- /dev/null +++ b/rts/gl/rtsSourceCode.h @@ -0,0 +1,65 @@ +#ifndef RTSSOURCECODE_H +#define RTSSOURCECODE_H + +#include +#include +#include +#include + +using namespace std; + +///This class defines generic source code that can be loaded from text files. It is primarily used by the rts_glShaderProgram class for GLSL programming. + +class rtsSourceCode +{ +public: + vector source; //the actual source code + void clear() /// +//#include "windows.h" +#include +#include "rtsSourceCode.h" + +class rts_glShaderObject +{ +private: + void init() + { + id = 0; + compiled = false; + type = GL_FRAGMENT_SHADER; + } +public: + bool compiled; + GLenum type; + rtsSourceCode source; + GLuint id; + string log; + + rts_glShaderObject(GLenum type, const char* filename) + { + init(); //initialize the shader + SetType(type); //set the shader type + LoadSource(filename); //load the source code + } + rts_glShaderObject(GLenum type, rtsSourceCode sourceCode) + { + init(); //initialize the shader + SetType(type); //set the shader type + source = sourceCode; + } + rts_glShaderObject() + { + init(); + } + rts_glShaderObject(GLenum type) + { + init(); + SetType(type); + } + void LoadSource(const char* filename) + { + source = rtsSourceCode(filename); //get the shader source code + + } + void SetType(GLenum type) + { + if(id != 0) //if a shader currently exists, delete it + { + glDeleteShader(id); + id = 0; + } + type = type; + id = glCreateShader(type); //create a shader object + if(id == 0) //if a shader was not created, log an error + { + log = "Error getting shader ID from OpenGL"; + return; + } + } + void UploadSource() + { + //create the structure for the shader source code + GLsizei count = source.source.size(); + GLchar** code_string = new GLchar*[count]; + GLint* length = new GLint[count]; + for(int l = 0; l + +using namespace std; + +class rts_glShaderProgram +{ +private: + void get_uniforms() + { + GLint num_uniforms; + glGetProgramiv(id, GL_ACTIVE_UNIFORMS, &num_uniforms); //get the number of uniform variables + GLint max_name_length; + glGetProgramiv(id, GL_ACTIVE_UNIFORM_MAX_LENGTH, &max_name_length); //get the maximum uniform name length + GLchar* name_buffer = new GLchar[max_name_length]; //create a buffer to store the name + GLsizei length; //I'm not using these yet + GLint size; + GLenum type; //variable's data type + GLint location; //GPU location of the variable + for(int i=0; i shader_list; //list of opengl shaders + vector uniform_list; //list of active uniform variables + vector texture_list; //list of texture maps + + rts_glShaderProgram() + { + linked = false; + id = 0; + } + void AttachShader(rts_glShaderObject shader) + { + if(id == 0) + { + Init(); + } + if(shader.id == 0) //if the shader is invalid + { + log = "Shader is invalid"; + return; + } + + //attach the shader to the program + glAttachShader(id, shader.id); //attach the shader to the program in OpenGL + CHECK_OPENGL_ERROR + shader_list.push_back(shader); //push the shader onto our list for later access + } + //type = GL_FRAGMENT_SHADER or GL_VERTEX_SHADER + void AttachShader(GLenum type, const char* filename) + { + rts_glShaderObject shader(type, filename); + AttachShader(shader); + } + void AttachShader(GLenum type, rtsSourceCode source) + { + rts_glShaderObject shader(type, source); + AttachShader(shader); + } + void PrintLog() + { + cout<::iterator iter; + for(iter = shader_list.begin(); iter != shader_list.end(); iter++) + { + (*iter).Compile(); + //(*iter).PrintLog(); + } + } + void Link() + { + glLinkProgram(id); //link the current shader program + GLint link_status; //test to see if the link went alright + glGetProgramiv(id, GL_LINK_STATUS, &link_status); + if(link_status != GL_TRUE) + { + linked = false; + } + else + linked = true; + + GLsizei length; + GLchar buffer[1000]; + glGetProgramInfoLog(id, 1000, &length, buffer); + log = buffer; + + get_uniforms(); //create the list of active uniform variables + } + void BeginProgram() + { + CHECK_OPENGL_ERROR + if(id == 0) //if the program is invalid, return + { + log = "Invalid program, cannot use."; + return; + } + if(!linked) + { + cout<<"Shader Program used without being linked."< 0) + glActiveTexture(GL_TEXTURE0); + CHECK_OPENGL_ERROR + + //return to OpenGL default shading + glUseProgram(0); + CHECK_OPENGL_ERROR + } + void PrintUniforms() + { + cout<<"Shader Uniforms: "< +#include + +using namespace std; + +enum rtsUniformEnum {RTS_FLOAT, RTS_INT, RTS_BOOL, RTS_FLOAT_MATRIX}; + +///This class stores a single uniform variable for GLSL and is designed to be used by the rts_glShaderProgram class. +struct rts_glShaderUniform +{ +public: + string name; //the name of the variable + GLint location; //the location in the program + void* p_value; //pointer to the global data representing the value in main memory + GLenum type; //variable type (float, int, vec2, etc.) + //rtsUniformEnum rts_type; //type of variable in rts format + //unsigned int num; //the number of values required by the variable (1 for float, 2 for vec2, etc.) + string log; + + //void convert_type(GLenum gl_type); //converts the OpenGL data type to something useful for rts + void submit_to_gpu() + { + if(location < 0) + return; + if(p_value == NULL) + { + cout<<"Error in uniform address: "< +#include "rts/math/vector.h" +#include "rts/gl/error.h" +#include + +namespace rts{ + +///This class stores an OpenGL texture map and is used by rts_glShaderProgram. +class glTexture +{ +private: + void get_type() //guesses the texture type based on the size + { + if(size[1] == 0) + texture_type = GL_TEXTURE_1D; + else if(size[2] == 0) + texture_type = GL_TEXTURE_2D; + else + texture_type = GL_TEXTURE_3D; + } + void set_wrapping() //set the texture wrapping based on the dimensions + { + CHECK_OPENGL_ERROR + switch(texture_type) + { + case GL_TEXTURE_3D: + glTexParameteri(texture_type, GL_TEXTURE_WRAP_R_EXT, GL_REPEAT); + case GL_TEXTURE_2D: + glTexParameteri(texture_type, GL_TEXTURE_WRAP_T, GL_MIRRORED_REPEAT); + case GL_TEXTURE_1D: + glTexParameteri(texture_type, GL_TEXTURE_WRAP_S, GL_REPEAT); + break; + case GL_TEXTURE_RECTANGLE_ARB: + glTexParameteri(texture_type, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(texture_type, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + break; + + default: + break; + } + CHECK_OPENGL_ERROR + } + //void set_bits(GLvoid* bits); +public: + vector size; //vector representing the size of the texture + GLuint name; //texture name assigned by OpenGL + GLenum texture_type; //1D, 2D, 3D + GLint internal_format; //number of components (ex. 4 for RGBA) + GLenum pixel_format; //type of data (RGBA, LUMINANCE) + GLenum data_type; //data type of the bits (float, int, etc.) + + //constructor + glTexture() + { + name = 0; + } + glTexture(GLvoid *bits, + GLenum type = GL_TEXTURE_2D, + GLsizei width = 256, + GLsizei height = 256, + GLsizei depth = 0, + GLint internalformat = 1, + GLenum format = GL_LUMINANCE, + GLenum datatype = GL_UNSIGNED_BYTE, + GLint interpolation = GL_LINEAR) + { + init(bits, type, width, height, depth, internalformat, format, datatype, interpolation); + } + + void begin() + { + glEnable(texture_type); + CHECK_OPENGL_ERROR + glBindTexture(texture_type, name); + CHECK_OPENGL_ERROR + } + void end() + { + glDisable(texture_type); + CHECK_OPENGL_ERROR + } + + ///Creates an OpenGL texture map. This function requires basic information about the texture map as well as a pointer to the bit data describing the texture. + void init(GLvoid *bits, + GLenum type = GL_TEXTURE_2D, + GLsizei width = 256, + GLsizei height = 256, + GLsizei depth = 0, + GLint internalformat = 1, + GLenum format = GL_LUMINANCE, + GLenum datatype = GL_UNSIGNED_BYTE, + GLint interpolation = GL_LINEAR) + { + CHECK_OPENGL_ERROR + if(name != 0) + glDeleteTextures(1, &name); + + + CHECK_OPENGL_ERROR + if(datatype == GL_FLOAT) + { + glPixelStorei(GL_PACK_ALIGNMENT, 4); + glPixelStorei(GL_UNPACK_ALIGNMENT, 4); //I honestly don't know what this does but it fixes problems + } + else if(datatype == GL_UNSIGNED_BYTE) + { + //glPixelStorei(GL_UNPACK_ALIGNMENT, 1); + //glPixelStorei(GL_PACK_ALIGNMENT, 1); + } + else if(datatype == GL_UNSIGNED_SHORT) + { + //glPixelStorei(GL_UNPACK_ALIGNMENT, 2); + //glPixelStorei(GL_PACK_ALIGNMENT, 2); + } + CHECK_OPENGL_ERROR + glGenTextures(1, &name); //get the texture name from OpenGL + //cout<<"OpenGL Name: "<(width, height, depth); //assign the texture size + //get_type(); //guess the type based on the size + texture_type = type; //set the type of texture + glEnable(texture_type); //enable the texture map + CHECK_OPENGL_ERROR + glBindTexture(texture_type, name); //bind the texture for editing + CHECK_OPENGL_ERROR + set_wrapping(); //set the texture wrapping parameters + CHECK_OPENGL_ERROR + glTexParameteri(texture_type, GL_TEXTURE_MAG_FILTER, interpolation); //set filtering + CHECK_OPENGL_ERROR + glTexParameteri(texture_type, GL_TEXTURE_MIN_FILTER, interpolation); + CHECK_OPENGL_ERROR + internal_format = internalformat; //set the number of components per pixel + pixel_format = format; //set the pixel format + data_type = datatype; //set the data type + SetBits(bits); //send the bits to the OpenGL driver + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //replace the specified vertex color + CHECK_OPENGL_ERROR + glDisable(texture_type); + + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE); + glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_EDGE); + } + void Clean() + { + if(name != 0) + { + glDeleteTextures(1, &name); + CHECK_OPENGL_ERROR + name = 0; + } + } + void SetBits(GLvoid *bits) + { + glEnable(texture_type); //enable the texture map + CHECK_OPENGL_ERROR + glBindTexture(texture_type, name); + CHECK_OPENGL_ERROR + + switch(texture_type) + { + case GL_TEXTURE_3D: + glTexImage3D(texture_type, 0, internal_format, size[0], size[1], size[2], 0, pixel_format, data_type, bits); + CHECK_OPENGL_ERROR + break; + case GL_TEXTURE_2D: + case GL_TEXTURE_RECTANGLE_ARB: + glTexImage2D(texture_type, 0, internal_format, size[0], size[1], 0, pixel_format, data_type, bits); + CHECK_OPENGL_ERROR + break; + case GL_TEXTURE_1D: + glTexImage1D(texture_type, 0, internal_format, size[0], 0, pixel_format, data_type, bits); + CHECK_OPENGL_ERROR + break; + default: + //glTexImage2D(texture_type, 0, internal_format, size.x, size.y, 0, pixel_format, data_type, bits); + break; + } + CHECK_OPENGL_ERROR + } + void ResetBits(GLvoid *bits) + { + glEnable(texture_type); //enable the texture map + CHECK_OPENGL_ERROR + glBindTexture(texture_type, name); + CHECK_OPENGL_ERROR + + switch(texture_type) + { + case GL_TEXTURE_3D: + //glTexImage3D(texture_type, 0, internal_format, size.x, size.y, size.z, 0, pixel_format, data_type, bits); + break; + case GL_TEXTURE_2D: + case GL_TEXTURE_RECTANGLE_ARB: + glTexSubImage2D(texture_type, 0, 0, 0, size[0], size[1], pixel_format, data_type, bits); + CHECK_OPENGL_ERROR + break; + case GL_TEXTURE_1D: + //glTexImage1D(texture_type, 0, internal_format, size.x, 0, pixel_format, data_type, bits); + break; + default: + //glTexImage2D(texture_type, 0, internal_format, size.x, size.y, 0, pixel_format, data_type, bits); + break; + } + glDisable(texture_type); + CHECK_OPENGL_ERROR + } + void* GetBits(GLenum format, GLenum type) + { + //returns the texture data + + int components; + switch(format) + { + case GL_RED: + case GL_GREEN: + case GL_BLUE: + case GL_ALPHA: + case GL_LUMINANCE: + components = 1; + break; + case GL_LUMINANCE_ALPHA: + components = 2; + break; + case GL_RGB: + case GL_BGR: + components = 3; + break; + case GL_RGBA: + case GL_BGRA: + components = 4; + break; + } + + int type_size; + switch(type) + { + case GL_UNSIGNED_BYTE: + case GL_BYTE: + type_size = sizeof(char); + break; + case GL_UNSIGNED_SHORT: + case GL_SHORT: + type_size = sizeof(short); + break; + case GL_UNSIGNED_INT: + case GL_INT: + type_size = sizeof(int); + break; + case GL_FLOAT: + type_size = sizeof(float); + break; + } + + //allocate memory for the texture + void* result = malloc(components*type_size * size[0] * size[1]); + + begin(); + glGetTexImage(texture_type, 0, format, type, result); + + CHECK_OPENGL_ERROR + end(); + + + return result; + + } +}; + +} + +#define RTS_UNKNOWN 0 + +#endif diff --git a/rts/math/complex.h b/rts/math/complex.h index 1a338f6..6c306bd 100644 --- a/rts/math/complex.h +++ b/rts/math/complex.h @@ -15,12 +15,12 @@ namespace rts { template -struct rtsComplex +struct complex { T r, i; //default constructor - CUDA_CALLABLE rtsComplex() + CUDA_CALLABLE complex() { r = 0.0; i = 0.0; @@ -49,16 +49,16 @@ struct rtsComplex } //constructor when given real and imaginary values - CUDA_CALLABLE rtsComplex(T r, T i) + CUDA_CALLABLE complex(T r, T i) { this->r = r; this->i = i; } //return the current value multiplied by i - CUDA_CALLABLE rtsComplex imul() + CUDA_CALLABLE complex imul() { - rtsComplex result; + complex result; result.r = -i; result.i = r; @@ -68,70 +68,70 @@ struct rtsComplex //ARITHMETIC OPERATORS-------------------- //binary + operator (returns the result of adding two complex values) - CUDA_CALLABLE rtsComplex operator+ (const rtsComplex rhs) + CUDA_CALLABLE complex operator+ (const complex rhs) { - rtsComplex result; + complex result; result.r = r + rhs.r; result.i = i + rhs.i; return result; } - CUDA_CALLABLE rtsComplex operator+ (const T rhs) + CUDA_CALLABLE complex operator+ (const T rhs) { - rtsComplex result; + complex result; result.r = r + rhs; result.i = i; return result; } //binary - operator (returns the result of adding two complex values) - CUDA_CALLABLE rtsComplex operator- (const rtsComplex rhs) + CUDA_CALLABLE complex operator- (const complex rhs) { - rtsComplex result; + complex result; result.r = r - rhs.r; result.i = i - rhs.i; return result; } //binary - operator (returns the result of adding two complex values) - CUDA_CALLABLE rtsComplex operator- (const T rhs) + CUDA_CALLABLE complex operator- (const T rhs) { - rtsComplex result; + complex result; result.r = r - rhs; result.i = i; return result; } //binary MULTIPLICATION operators (returns the result of multiplying complex values) - CUDA_CALLABLE rtsComplex operator* (const rtsComplex rhs) + CUDA_CALLABLE complex operator* (const complex rhs) { - rtsComplex result; + complex result; result.r = r * rhs.r - i * rhs.i; result.i = r * rhs.i + i * rhs.r; return result; } - CUDA_CALLABLE rtsComplex operator* (const T rhs) + CUDA_CALLABLE complex operator* (const T rhs) { - return rtsComplex(r * rhs, i * rhs); + return complex(r * rhs, i * rhs); } //binary DIVISION operators (returns the result of dividing complex values) - CUDA_CALLABLE rtsComplex operator/ (const rtsComplex rhs) + CUDA_CALLABLE complex operator/ (const complex rhs) { - rtsComplex result; + complex result; T denom = rhs.r * rhs.r + rhs.i * rhs.i; result.r = (r * rhs.r + i * rhs.i) / denom; result.i = (- r * rhs.i + i * rhs.r) / denom; return result; } - CUDA_CALLABLE rtsComplex operator/ (const T rhs) + CUDA_CALLABLE complex operator/ (const T rhs) { - return rtsComplex(r / rhs, i / rhs); + return complex(r / rhs, i / rhs); } //ASSIGNMENT operators----------------------------------- - CUDA_CALLABLE rtsComplex & operator=(const rtsComplex &rhs) + CUDA_CALLABLE complex & operator=(const complex &rhs) { //check for self-assignment if(this != &rhs) @@ -141,7 +141,7 @@ struct rtsComplex } return *this; } - CUDA_CALLABLE rtsComplex & operator=(const T &rhs) + CUDA_CALLABLE complex & operator=(const T &rhs) { this->r = rhs; this->i = 0; @@ -150,34 +150,34 @@ struct rtsComplex } //arithmetic assignment operators - CUDA_CALLABLE rtsComplex operator+=(const rtsComplex &rhs) + CUDA_CALLABLE complex operator+=(const complex &rhs) { *this = *this + rhs; return *this; } - CUDA_CALLABLE rtsComplex operator+=(const T &rhs) + CUDA_CALLABLE complex operator+=(const T &rhs) { *this = *this + rhs; return *this; } - CUDA_CALLABLE rtsComplex operator*=(const rtsComplex &rhs) + CUDA_CALLABLE complex operator*=(const complex &rhs) { *this = *this * rhs; return *this; } - CUDA_CALLABLE rtsComplex operator*=(const T &rhs) + CUDA_CALLABLE complex operator*=(const T &rhs) { *this = *this * rhs; return *this; } //divide and assign - CUDA_CALLABLE rtsComplex operator/=(const rtsComplex &rhs) + CUDA_CALLABLE complex operator/=(const complex &rhs) { *this = *this / rhs; return *this; } - CUDA_CALLABLE rtsComplex operator/=(const T &rhs) + CUDA_CALLABLE complex operator/=(const T &rhs) { *this = *this / rhs; return *this; @@ -189,9 +189,9 @@ struct rtsComplex return std::sqrt(r * r + i * i); } - CUDA_CALLABLE rtsComplex log() + CUDA_CALLABLE complex log() { - rtsComplex result; + complex result; result.r = std::log(std::sqrt(r * r + i * i)); result.i = std::atan2(i, r); @@ -199,9 +199,9 @@ struct rtsComplex return result; } - CUDA_CALLABLE rtsComplex exp() + CUDA_CALLABLE complex exp() { - rtsComplex result; + complex result; T e_r = std::exp(r); result.r = e_r * std::cos(i); @@ -216,18 +216,18 @@ struct rtsComplex return pow((double)y); }*/ - CUDA_CALLABLE rtsComplex pow(T y) + CUDA_CALLABLE complex pow(T y) { - rtsComplex result; + complex result; result = log() * y; return result.exp(); } - CUDA_CALLABLE rtsComplex sqrt() + CUDA_CALLABLE complex sqrt() { - rtsComplex result; + complex result; //convert to polar coordinates T a = std::sqrt(r*r + i*i); @@ -253,7 +253,7 @@ struct rtsComplex } //COMPARISON operators - CUDA_CALLABLE bool operator==(rtsComplex rhs) + CUDA_CALLABLE bool operator==(complex rhs) { if(r == rhs.r && i == rhs.i) return true; @@ -267,72 +267,44 @@ struct rtsComplex return false; } - /*//FRIEND functions - //unary minus operator (for negating the complex number) - template CUDA_CALLABLE friend complex operator-(const complex &rhs); - - //multiplication by T values when the complex number isn't on the left hand side - template CUDA_CALLABLE friend complex operator*(const A a, const complex b); - - //division by T values when the complex number isn't on the left hand side - template CUDA_CALLABLE friend complex operator/(const A a, const complex b); - - //POW function - //template CUDA_CALLABLE friend complex pow(const complex x, T y); - template CUDA_CALLABLE friend complex pow(const complex x, int y); - - //log function - template CUDA_CALLABLE friend complex log(complex x); - - //exp function - template CUDA_CALLABLE friend complex exp(complex x); - - //sqrt function - template CUDA_CALLABLE friend complex sqrt(complex x); - - //trigonometric functions - template CUDA_CALLABLE friend complex sin(complex x); - - template CUDA_CALLABLE friend complex cos(complex x);*/ - }; } //end RTS namespace //addition template -CUDA_CALLABLE static rts::rtsComplex operator+(const double a, const rts::rtsComplex b) +CUDA_CALLABLE static rts::complex operator+(const double a, const rts::complex b) { - return rts::rtsComplex(a + b.r, b.i); + return rts::complex(a + b.r, b.i); } //subtraction with a real value template -CUDA_CALLABLE static rts::rtsComplex operator-(const double a, const rts::rtsComplex b) +CUDA_CALLABLE static rts::complex operator-(const double a, const rts::complex b) { - return rts::rtsComplex(a - b.r, -b.i); + return rts::complex(a - b.r, -b.i); } //minus sign template -CUDA_CALLABLE static rts::rtsComplex operator-(const rts::rtsComplex &rhs) +CUDA_CALLABLE static rts::complex operator-(const rts::complex &rhs) { - return rts::rtsComplex(-rhs.r, -rhs.i); + return rts::complex(-rhs.r, -rhs.i); } //multiply a T value by a complex value template -CUDA_CALLABLE static rts::rtsComplex operator*(const double a, const rts::rtsComplex b) +CUDA_CALLABLE static rts::complex operator*(const double a, const rts::complex b) { - return rts::rtsComplex(a * b.r, a * b.i); + return rts::complex((T)a * b.r, (T)a * b.i); } //divide a T value by a complex value template -CUDA_CALLABLE static rts::rtsComplex operator/(const double a, const rts::rtsComplex b) +CUDA_CALLABLE static rts::complex operator/(const double a, const rts::complex b) { //return complex(a * b.r, a * b.i); - rts::rtsComplex result; + rts::complex result; T denom = b.r * b.r + b.i * b.i; @@ -350,41 +322,41 @@ CUDA_CALLABLE static complex pow(complex x, int y) }*/ template -CUDA_CALLABLE static rts::rtsComplex pow(rts::rtsComplex x, T y) +CUDA_CALLABLE static rts::complex pow(rts::complex x, T y) { return x.pow(y); } //log function template -CUDA_CALLABLE static rts::rtsComplex log(rts::rtsComplex x) +CUDA_CALLABLE static rts::complex log(rts::complex x) { return x.log(); } //exp function template -CUDA_CALLABLE static rts::rtsComplex exp(rts::rtsComplex x) +CUDA_CALLABLE static rts::complex exp(rts::complex x) { return x.exp(); } //sqrt function template -CUDA_CALLABLE static rts::rtsComplex sqrt(rts::rtsComplex x) +CUDA_CALLABLE static rts::complex sqrt(rts::complex x) { return x.sqrt(); } template -CUDA_CALLABLE static T abs(rts::rtsComplex a) +CUDA_CALLABLE static T abs(rts::complex a) { return a.abs(); } template -CUDA_CALLABLE static T real(rts::rtsComplex a) +CUDA_CALLABLE static T real(rts::complex a) { return a.r; } @@ -396,16 +368,16 @@ CUDA_CALLABLE static float real(float a) } template -CUDA_CALLABLE static T imag(rts::rtsComplex a) +CUDA_CALLABLE static T imag(rts::complex a) { return a.i; } //trigonometric functions template -CUDA_CALLABLE rts::rtsComplex sin(const rts::rtsComplex x) +CUDA_CALLABLE rts::complex sin(const rts::complex x) { - rts::rtsComplex result; + rts::complex result; result.r = std::sin(x.r) * std::cosh(x.i); result.i = std::cos(x.r) * std::sinh(x.i); @@ -413,9 +385,9 @@ CUDA_CALLABLE rts::rtsComplex sin(const rts::rtsComplex x) } template -CUDA_CALLABLE rts::rtsComplex cos(const rts::rtsComplex x) +CUDA_CALLABLE rts::complex cos(const rts::complex x) { - rts::rtsComplex result; + rts::complex result; result.r = std::cos(x.r) * std::cosh(x.i); result.i = -(std::sin(x.r) * std::sinh(x.i)); @@ -424,12 +396,16 @@ CUDA_CALLABLE rts::rtsComplex cos(const rts::rtsComplex x) template -std::ostream& operator<<(std::ostream& os, rts::rtsComplex x) +std::ostream& operator<<(std::ostream& os, rts::complex x) { os< 3 && __GNUC_MINOR__ > 7 +template using rtsComplex = rts::complex; +#endif + #endif diff --git a/rts/math/function.h b/rts/math/function.h new file mode 100644 index 0000000..d3c200e --- /dev/null +++ b/rts/math/function.h @@ -0,0 +1,123 @@ +#ifndef RTS_FUNCTION_H +#define RTS_FUNCTION_H + +#include + +namespace rts{ + +template +class function +{ + //datapoint class for storing function points + struct dataPoint + { + X x; + Y y; + }; + + //function data + std::vector f; + + //comparison function for searching lambda + static bool findCeiling(dataPoint a, dataPoint b) + { + return (a.x > b.x); + } + + +public: + Y linear(X x) + { + //declare an iterator + typename std::vector< dataPoint >::iterator it; + + dataPoint s; + s.x = x; + + it = search(f.begin(), f.end(), &s, &s + 1, &function::findCeiling); + + //if the wavelength is past the end of the list, return the back + if(it == f.end()) + return f.back().y; + //if the wavelength is before the beginning of the list, return the front + else if(it == f.begin()) + return f.front().y; + //otherwise interpolate + else + { + X xMax = (*it).x; + X xMin = (*(it - 1)).x; + //std::cout<::iterator it; + + dataPoint s; + s.x = x; + s.y = y; + + it = search(f.begin(), f.end(), &s, &s + 1, &function::findCeiling); + + //if the function value is past the end of the vector, add it to the back + if(it == f.end()) + return f.push_back(s); + //otherwise add the value at the iterator position + else + { + f.insert(it, s); + } + + } + + X getX(unsigned int i) + { + return f[i].x; + } + + Y getY(unsigned int i) + { + return f[i].y; + } + + unsigned int getN() + { + return f.size(); + } + + dataPoint operator[](int i) + { + return f[i]; + } + + function operator+(Y r) + { + function result; + + //add r to every point in f + for(int i=0; i -struct rtsMatrix +struct matrix { //the matrix will be stored in column-major order (compatible with OpenGL) T M[N*N]; - rtsMatrix() + matrix() { for(int r=0; r operator=(T rhs) + matrix operator=(T rhs) { int Nsq = N*N; for(int i=0; i operator=(rtsMatrix rhs) + /*matrix operator=(matrix rhs) { for(int i=0; i operator*(rtsVector rhs) + vector operator*(vector rhs) { - rtsVector result; + vector result; for(int r=0; r -std::ostream& operator<<(std::ostream& os, rts::rtsMatrix M) +std::ostream& operator<<(std::ostream& os, rts::matrix M) { os< 3 && __GNUC_MINOR__ > 7 +template using rtsMatrix = rts::matrix; +#endif + #endif diff --git a/rts/math/point.h b/rts/math/point.h index 801347f..bd9c4cb 100644 --- a/rts/math/point.h +++ b/rts/math/point.h @@ -1,7 +1,7 @@ #ifndef RTS_rtsPoint_H #define RTS_rtsPoint_H -#include "rts/math/vector.h" +#include "rts/math/vector.h" #include #include "rts/cuda/callable.h" @@ -9,17 +9,17 @@ namespace rts { template -struct rtsPoint +struct point { T p[N]; - CUDA_CALLABLE rtsPoint() - { + CUDA_CALLABLE point() + { } //efficiency constructor, makes construction easier for 1D-4D vectors - CUDA_CALLABLE rtsPoint(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0) + CUDA_CALLABLE point(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0) { if(N >= 1) p[0] = x; @@ -29,51 +29,51 @@ struct rtsPoint p[2] = z; if(N >= 4) p[3] = w; - } - - //arithmetic operators - CUDA_CALLABLE rts::rtsPoint operator+(rts::rtsVector v) - { - rts::rtsPoint r; - - //calculate the position of the resulting rtsPoint - for(int i=0; i operator-(rts::rtsVector v) - { - rts::rtsPoint r; - - //calculate the position of the resulting rtsPoint - for(int i=0; i operator-(rts::rtsPoint rhs) - { - rts::rtsVector r; - - //calculate the position of the resulting rtsPoint - for(int i=0; i operator*(T rhs) - { - rts::rtsPoint r; - - //calculate the position of the resulting rtsPoint - for(int i=0; i operator+(vector v) + { + rts::point r; + + //calculate the position of the resulting point + for(int i=0; i operator-(vector v) + { + rts::point r; + + //calculate the position of the resulting point + for(int i=0; i operator-(point rhs) + { + vector r; + + //calculate the position of the resulting point + for(int i=0; i operator*(T rhs) + { + rts::point r; + + //calculate the position of the resulting point + for(int i=0; i -std::ostream& operator<<(std::ostream& os, rts::rtsPoint p) -{ - os< -CUDA_CALLABLE rts::rtsPoint operator*(T lhs, rts::rtsPoint rhs) -{ - rts::rtsPoint r; - - return rhs * lhs; -} - - - -#endif +}; + +} //end namespace rts + +template +std::ostream& operator<<(std::ostream& os, rts::point p) +{ + os< +CUDA_CALLABLE rts::point operator*(T lhs, rts::point rhs) +{ + rts::point r; + + return rhs * lhs; +} + +#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 +template using rtsPoint = rts::point; +#endif + +#endif diff --git a/rts/math/quad.h b/rts/math/quad.h index cd1ccd3..3913358 100644 --- a/rts/math/quad.h +++ b/rts/math/quad.h @@ -6,13 +6,14 @@ #include "rts/math/vector.h" #include "rts/math/point.h" #include "rts/math/triangle.h" +#include "rts/math/quaternion.h" #include namespace rts{ -//template for a rtsQuadangle class in ND space +//template for a quadangle class in ND space template -struct rtsQuad +struct quad { /* C------------------>O @@ -28,17 +29,17 @@ struct rtsQuad T B[N]; T C[N];*/ - rts::rtsPoint A; - rts::rtsVector X; - rts::rtsVector Y; + rts::point A; + rts::vector X; + rts::vector Y; - CUDA_CALLABLE rtsQuad() + CUDA_CALLABLE quad() { } - CUDA_CALLABLE rtsQuad(rtsPoint a, rtsPoint b, rtsPoint c) + CUDA_CALLABLE quad(point a, point b, point c) { A = a; @@ -47,48 +48,90 @@ struct rtsQuad } - CUDA_CALLABLE rtsQuad(rts::rtsPoint pMin, rts::rtsPoint pMax, rts::rtsVector normal) + /**************************************************************** + Constructor - create a quad from two points and a normal + ****************************************************************/ + CUDA_CALLABLE quad(rts::point pMin, rts::point pMax, rts::vector normal) { - //assign the corner rtsPoint + //assign the corner point A = pMin; //compute the vector from pMin to pMax - rts::rtsVector v0; + rts::vector v0; v0 = pMax - pMin; //compute the cross product of A and the plane normal - rts::rtsVector v1; + rts::vector v1; v1 = v0.cross(normal); - //calculate rtsPoint B - rts::rtsPoint B; + //calculate point B + rts::point B; B = A + v0 * 0.5 + v1 * 0.5; //calculate rtsPoint C - rts::rtsPoint C; + rts::point C; C = A + v0 * 0.5 - v1 * 0.5; //calculate X and Y X = B - A; Y = C - A; + } + + /******************************************************************* + Constructor - create a quad from a position, normal, and rotation + *******************************************************************/ + CUDA_CALLABLE quad(rts::point c, rts::vector normal, T width, T height, T theta) + { + + //compute the X direction - start along world-space X + Y = rts::vector(0, 1, 0); + if(Y == normal) + Y = rts::vector(0, 0, 1); + + X = Y.cross(normal).norm(); + + std::cout< q; + q.CreateRotation(theta, normal); + X = q.toMatrix3() * X; + Y = normal.cross(X); + //normalize everything + X = X.norm(); + Y = Y.norm(); + //scale to match the quad width and height + X = X * width; + Y = Y * height; + //set the corner of the plane + A = c - X * 0.5 - Y * 0.5; + std::cout< p(T a, T b) + /******************************************* + Return the normal for the quad + *******************************************/ + CUDA_CALLABLE rts::vector n() { - rts::rtsPoint result; + return (X.cross(Y)).norm(); + } + + CUDA_CALLABLE rts::point p(T a, T b) + { + rts::point result; //given the two parameters a, b = [0 1], returns the position in world space result = A + X * a + Y * b; return result; } - CUDA_CALLABLE rts::rtsPoint operator()(T a, T b) + CUDA_CALLABLE rts::point operator()(T a, T b) { return p(a, b); } @@ -106,15 +149,15 @@ struct rtsQuad } - CUDA_CALLABLE rtsQuad operator*(T rhs) + CUDA_CALLABLE quad operator*(T rhs) { //scales the plane by a scalar value - //compute the center rtsPoint - rts::rtsPoint c = A + X*0.5 + Y*0.5; + //compute the center point + rts::point c = A + X*0.5 + Y*0.5; - //create the new rtsQuadangle - rtsQuad result; + //create the new quadangle + quad result; result.X = X * rhs; result.Y = Y * rhs; result.A = c - result.X*0.5 - result.Y*0.5; @@ -123,13 +166,13 @@ struct rtsQuad } - CUDA_CALLABLE T dist(rtsPoint p) + CUDA_CALLABLE T dist(point p) { //compute the distance between a point and this quad //first break the quad up into two triangles - rtsTriangle T0(A, A+X, A+Y); - rtsTriangle T1(A+X+Y, A+X, A+Y); + triangle T0(A, A+X, A+Y); + triangle T1(A+X+Y, A+X, A+Y); ptype d0 = T0.dist(p); @@ -140,12 +183,22 @@ struct rtsQuad else return d1; } + + CUDA_CALLABLE T dist_max(point p) + { + T da = (A - p).len(); + T db = (A+X - p).len(); + T dc = (A+Y - p).len(); + T dd = (A+X+Y - p).len(); + + return max( da, max(db, max(dc, dd) ) ); + } }; } //end namespace rts template -std::ostream& operator<<(std::ostream& os, rts::rtsQuad R) +std::ostream& operator<<(std::ostream& os, rts::quad R) { os< -class rtsQuaternion +class quaternion { public: T w; @@ -16,18 +16,19 @@ public: void normalize(); void CreateRotation(T theta, T axis_x, T axis_y, T axis_z); - void CreateRotation(T theta, rtsVector axis); - rtsQuaternion operator*(rtsQuaternion &rhs); - rtsMatrix toMatrix(); + void CreateRotation(T theta, vector axis); + quaternion operator*(quaternion &rhs); + matrix toMatrix3(); + matrix toMatrix4(); - rtsQuaternion(); - rtsQuaternion(T w, T x, T y, T z); + quaternion(); + quaternion(T w, T x, T y, T z); }; template -void rtsQuaternion::normalize() +void quaternion::normalize() { double length=sqrt(w*w + x*x + y*y + z*z); w=w/length; @@ -37,23 +38,23 @@ void rtsQuaternion::normalize() } template -void rtsQuaternion::CreateRotation(T theta, T axis_x, T axis_y, T axis_z) +void quaternion::CreateRotation(T theta, T axis_x, T axis_y, T axis_z) { //assign the given Euler rotation to this quaternion - w = cos(theta/2.0); - x = axis_x*sin(theta/2.0); - y = axis_y*sin(theta/2.0); - z = axis_z*sin(theta/2.0); + w = (T)cos(theta/2.0); + x = axis_x*(T)sin(theta/2.0); + y = axis_y*(T)sin(theta/2.0); + z = axis_z*(T)sin(theta/2.0); } template -void rtsQuaternion::CreateRotation(T theta, rtsVector axis) +void quaternion::CreateRotation(T theta, vector axis) { CreateRotation(theta, axis[0], axis[1], axis[2]); } template -rtsQuaternion rtsQuaternion::operator *(rtsQuaternion ¶m) +quaternion quaternion::operator *(quaternion ¶m) { float A, B, C, D, E, F, G, H; @@ -67,7 +68,7 @@ rtsQuaternion rtsQuaternion::operator *(rtsQuaternion ¶m) G = (w + y)*(param.w - param.z); H = (w - y)*(param.w + param.z); - rtsQuaternion result; + quaternion result; result.w = B + (-E - F + G + H) /2; result.x = A - (E + F + G + H)/2; result.y = C + (E - F + G - H)/2; @@ -77,12 +78,12 @@ rtsQuaternion rtsQuaternion::operator *(rtsQuaternion ¶m) } template -rtsMatrix rtsQuaternion::toMatrix() +matrix quaternion::toMatrix3() { - rtsMatrix result; + matrix result; - double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2; + T wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2; // calculate coefficients @@ -92,62 +93,76 @@ rtsMatrix rtsQuaternion::toMatrix() yy = y * y2; yz = y * z2; zz = z * z2; wx = w * x2; wy = w * y2; wz = w * z2; - result(0, 0) = 1.0 - (yy + zz); + result(0, 0) = (T)1.0 - (yy + zz); result(0, 1) = xy - wz; - //m[0][0] = 1.0 - (yy + zz); m[1][0] = xy - wz; result(0, 2) = xz + wy; - //result(0, 3) = 0.0; - //m[2][0] = xz + wy; m[3][0] = 0.0; result(1, 0) = xy + wz; - result(1, 1) = 1.0 - (xx + zz); - //m[0][1] = xy + wz; m[1][1] = 1.0 - (xx + zz); + result(1, 1) = (T)1.0 - (xx + zz); result(1, 2) = yz - wx; - //result(1, 3) = 0.0; - //m[2][1] = yz - wx; m[3][1] = 0.0; result(2, 0) = xz - wy; result(2, 1) = yz + wx; - //m[0][2] = xz - wy; m[1][2] = yz + wx; - result(2, 2) = 1.0 - (xx + yy); - //result(3, 2) = 0.0; - //m[2][2] = 1.0 - (xx + yy); m[3][2] = 0.0; + result(2, 2) = (T)1.0 - (xx + yy); + + return result; +} + +template +matrix quaternion::toMatrix4() +{ + matrix result; + + + T wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2; - /*result(3, 0) = 0.0; - result(3, 1) = 0.0; - result(3, 2) = 0.0; - result(3, 3) = 1.0;*/ - //m[0][3] = 0; m[1][3] = 0; - //m[2][3] = 0; m[3][3] = 1; - /* - double* orientationmatrix=(double*)m; - char c; + // calculate coefficients + x2 = x + x; y2 = y + y; + z2 = z + z; + xx = x * x2; xy = x * y2; xz = x * z2; + yy = y * y2; yz = y * z2; zz = z * z2; + wx = w * x2; wy = w * y2; wz = w * z2; + + result(0, 0) = (T)1.0 - (yy + zz); + result(0, 1) = xy - wz; + + result(0, 2) = xz + wy; + + result(1, 0) = xy + wz; + result(1, 1) = (T)1.0 - (xx + zz); + + result(1, 2) = yz - wx; + + result(2, 0) = xz - wy; + result(2, 1) = yz + wx; - double* result=new double[16]; - double* array=(double*)m; - for(int i=0; i<16; i++) - result[i]=array[i]; - */ + result(2, 2) = (T)1.0 - (xx + yy); + + result(3, 3) = (T)1.0; return result; } template -rtsQuaternion::rtsQuaternion() +quaternion::quaternion() { w=0.0; x=0.0; y=0.0; z=0.0; } template -rtsQuaternion::rtsQuaternion(T c, T i, T j, T k) +quaternion::quaternion(T c, T i, T j, T k) { w=c; x=i; y=j; z=k; } -} //end rts namespace +} //end rts namespace + +#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 +template using rtsQuaternion = rts::quaternion; +#endif #endif diff --git a/rts/math/spherical_bessel.h b/rts/math/spherical_bessel.h index e846e0b..c0112be 100644 --- a/rts/math/spherical_bessel.h +++ b/rts/math/spherical_bessel.h @@ -10,7 +10,7 @@ namespace rts{ #define RTS_BESSEL_MAXIMUM_FLOAT -1e33 template -CUDA_CALLABLE void sbesselj(int n, rtsComplex x, rtsComplex* j) +CUDA_CALLABLE void sbesselj(int n, complex x, complex* j) { //compute the first bessel function if(n >= 0) @@ -36,7 +36,7 @@ CUDA_CALLABLE void sbesselj(int n, rtsComplex x, rtsComplex* j) } template -CUDA_CALLABLE void sbessely(int n, rtsComplex x, rtsComplex* y) +CUDA_CALLABLE void sbessely(int n, complex x, complex* y) { //compute the first bessel function if(n >= 0) @@ -56,14 +56,14 @@ CUDA_CALLABLE void sbessely(int n, rtsComplex x, rtsComplex* y) //spherical Hankel functions of the first kind template -CUDA_CALLABLE void sbesselh1(int n, rtsComplex x, rtsComplex* h) +CUDA_CALLABLE void sbesselh1(int n, complex x, complex* h) { //compute j_0 and j_1 - rtsComplex j[2]; + complex j[2]; sbesselj(1, x, j); //compute y_0 and y_1 - rtsComplex y[2]; + complex y[2]; sbessely(1, x, y); //compute the first-order Hhankel function diff --git a/rts/math/triangle.h b/rts/math/triangle.h new file mode 100644 index 0000000..1db7466 --- /dev/null +++ b/rts/math/triangle.h @@ -0,0 +1,189 @@ +#ifndef RTS_TRIANGLE_H +#define RTS_TRIANGLE_H + +//enable CUDA_CALLABLE macro +#include "rts/cuda/callable.h" +#include "rts/math/vector.h" +#include "rts/math/point.h" +#include + +namespace rts{ + +template +struct triangle +{ + /* + A------>B + | / + | / + | / + | / + | / + | / + C + */ + private: + + point A; + point B; + point C; + + CUDA_CALLABLE point _p(T s, T t) + { + //This function returns the point specified by p = A + s(B-A) + t(C-A) + vector E0 = B-A; + vector E1 = C-A; + + return A + s*E0 + t*E1; + } + + + public: + + + + CUDA_CALLABLE triangle() + { + + } + + CUDA_CALLABLE triangle(point a, point b, point c) + { + A = a; + B = b; + C = c; + } + + CUDA_CALLABLE rts::point operator()(T s, T t) + { + return _p(s, t); + } + + CUDA_CALLABLE point nearest(point p) + { + //comptue the distance between a point and this triangle + // This code is adapted from: http://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf + + vector E0 = B-A; + vector E1 = C-A; + vector D = A - p; + + T a = E0.dot(E0); + T b = E0.dot(E1); + T c = E1.dot(E1); + T d = E0.dot(D); + T e = E1.dot(D); + //T f = D.dot(D); + + T det = a*c - b*b; + T s = b*e - c*d; + T t = b*d - a*e; + + /*std::cout<<"E0: "<= 0 ? 0 : ( -e >= c ? 1 : -e/c ) ); + //done + } + } + else if(t < 0) + { + //region 5 + //std::cout<<"Region 5"<= 0 ? 0 : ( -d >= a ? 1 : -d/a ) ); + t = 0; + //done + } + else + { + //region 0 + //std::cout<<"Region 0"<= denom ? 1 : numer/denom ); + } + t = 1 - s; + //done + } + } + + //std::cout<<"s: "< p) + { + point n = nearest(p); + + return (p - n).len(); + } +}; + +} + +#endif diff --git a/rts/math/vector.h b/rts/math/vector.h index a9d77d2..e7a234b 100644 --- a/rts/math/vector.h +++ b/rts/math/vector.h @@ -1,8 +1,8 @@ #ifndef RTS_VECTOR_H #define RTS_VECTOR_H -#include -#include +#include +#include #include //#include "rts/point.h" #include "rts/cuda/callable.h" @@ -13,11 +13,11 @@ namespace rts template -struct rtsVector +struct vector { T v[N]; - CUDA_CALLABLE rtsVector() + CUDA_CALLABLE vector() { //memset(v, 0, sizeof(T) * N); for(int i=0; i= 1) v[0] = x; @@ -37,7 +37,7 @@ struct rtsVector v[3] = w; } - CUDA_CALLABLE rtsVector(const T(&data)[N]) + CUDA_CALLABLE vector(const T(&data)[N]) { memcpy(v, data, sizeof(T) * N); } @@ -54,25 +54,25 @@ struct rtsVector } - CUDA_CALLABLE rtsVector cart2sph() + CUDA_CALLABLE vector cart2sph() { //convert the vector from cartesian to spherical coordinates //x, y, z -> r, theta, phi (where theta = 0 to 2*pi) - rtsVector sph; - sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); - sph[1] = std::atan2(v[1], v[0]); + vector sph; + sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]); + sph[1] = std::atan2(v[1], v[0]); sph[2] = std::acos(v[2] / sph[0]); return sph; } - CUDA_CALLABLE rtsVector sph2cart() + CUDA_CALLABLE vector sph2cart() { //convert the vector from cartesian to spherical coordinates //r, theta, phi -> x, y, z (where theta = 0 to 2*pi) - rtsVector cart; + vector cart; cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]); cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]); cart[2] = v[0] * std::cos(v[2]); @@ -80,10 +80,10 @@ struct rtsVector return cart; } - CUDA_CALLABLE rtsVector norm() + CUDA_CALLABLE vector norm() { //compute and return the vector norm - rtsVector result; + vector result; //compute the vector length T l = len(); @@ -97,9 +97,9 @@ struct rtsVector return result; } - CUDA_CALLABLE rtsVector cross(rtsVector rhs) + CUDA_CALLABLE vector cross(vector rhs) { - rtsVector result; + vector result; //compute the cross product (only valid for 3D vectors) result[0] = v[1] * rhs[2] - v[2] * rhs[1]; @@ -109,7 +109,7 @@ struct rtsVector return result; } - CUDA_CALLABLE T dot(rtsVector rhs) + CUDA_CALLABLE T dot(vector rhs) { T result = (T)0; @@ -121,41 +121,49 @@ struct rtsVector } //arithmetic - CUDA_CALLABLE rtsVector operator+(rtsVector rhs) + CUDA_CALLABLE vector operator+(vector rhs) { - rtsVector result; + vector result; for(int i=0; i operator-(rtsVector rhs) + CUDA_CALLABLE vector operator-(vector rhs) { - rtsVector result; + vector result; for(int i=0; i operator*(T rhs) + CUDA_CALLABLE vector operator*(T rhs) { - rtsVector result; + vector result; for(int i=0; i operator/(T rhs) + CUDA_CALLABLE vector operator/(T rhs) { - rtsVector result; + vector result; for(int i=0; i rhs) + { + if ( (rhs.v[0] == v[0]) && (rhs.v[1] == v[1]) && (rhs.v[2] == v[2]) ) + return true; + + return false; } std::string toStr() @@ -186,7 +194,7 @@ struct rtsVector } //end namespace rts template -std::ostream& operator<<(std::ostream& os, rts::rtsVector v) +std::ostream& operator<<(std::ostream& os, rts::vector v) { os< v) //arithmetic operators template -CUDA_CALLABLE rts::rtsVector operator-(rts::rtsVector v) +CUDA_CALLABLE rts::vector operator-(rts::vector v) { - rts::rtsVector r; + rts::vector r; //negate the vector for(int i=0; i operator-(rts::rtsVector v) } template -CUDA_CALLABLE rts::rtsVector operator*(T lhs, rts::rtsVector rhs) +CUDA_CALLABLE rts::vector operator*(T lhs, rts::vector rhs) { - rts::rtsVector r; + rts::vector r; return rhs * lhs; } +#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 +template using rtsVector = rts::vector; +#endif + #endif diff --git a/rts/optics/material.h b/rts/optics/material.h index d2ff2f8..d6a4a99 100644 --- a/rts/optics/material.h +++ b/rts/optics/material.h @@ -9,6 +9,7 @@ #include #include #include "rts/math/complex.h" +#include "rts/math/function.h" #define PI 3.14159 @@ -49,7 +50,7 @@ namespace rts{ { //wavelength (in microns) T lambda; - rtsComplex n; + complex n; }; template @@ -111,10 +112,10 @@ namespace rts{ //read the entry from an input string - for(int i=0; i '9') { //cout<<"bad char"< > dispersion; + //std::vector< refIndex > dispersion; + function< T, complex > dispersion; //average refractive index (approximately 1.4) T n0; - void add(refIndex ri) + /*void add(refIndex ri) { //refIndex converted = convert(ri, measurement); dispersion.push_back(ri); - } - - void add(T lambda, rtsComplex n) - { - refIndex ri; - ri.lambda = lambda; - ri.n = n; - - dispersion.push_back(ri); - } + }*/ //comparison function for sorting static bool compare(refIndex a, refIndex b) @@ -222,12 +215,17 @@ namespace rts{ } //comparison function for searching lambda - static bool findCeiling(refIndex a, refIndex b) + /*static bool findCeiling(refIndex a, refIndex b) { return (a.lambda > b.lambda); - } + }*/ public: + void add(T lambda, complex n) + { + dispersion.insert(lambda, n); + } + std::string getName() { return name; @@ -244,6 +242,11 @@ namespace rts{ { n0 = n; } + + void setM(function< T, complex > m) + { + dispersion = m; + } unsigned int nSamples() { return dispersion.size(); @@ -268,9 +271,9 @@ namespace rts{ #ifdef FFTW_AVAILABLE //allocate memory for the FFT - rtsComplex* Chi2 = (rtsComplex*)fftw_malloc(sizeof(rtsComplex) * N * pf); - rtsComplex* Chi2FFT = (rtsComplex*)fftw_malloc(sizeof(rtsComplex) * N * pf); - rtsComplex* Chi1 = (rtsComplex*)fftw_malloc(sizeof(rtsComplex) * N * pf); + complex* Chi2 = (complex*)fftw_malloc(sizeof(complex) * N * pf); + complex* Chi2FFT = (complex*)fftw_malloc(sizeof(complex) * N * pf); + complex* Chi1 = (complex*)fftw_malloc(sizeof(complex) * N * pf); //create an FFT plan for the forward and inverse transforms fftw_plan planForward, planInverse; @@ -301,8 +304,8 @@ namespace rts{ } //use linear interpolation between the start and end points to pad the spectrum - //rtsComplex nMin = dispersion.back(); - //rtsComplex nMax = dispersion.front(); + //complex nMin = dispersion.back(); + //complex nMax = dispersion.front(); T a; for(int i=N; i j(0, 1); + complex j(0, 1); for(int i=0; i(0.0, k)); + /*//create a default refractive index refIndex def; def.lambda = lambda; def.n.real(n); def.n.imag(k); add(def); - + */ //set n0 n0 = n; } @@ -419,9 +423,9 @@ namespace rts{ void fromFile(std::string filename, std::string format = "", T scaleA = 1.0) { - name = filename; + name = filename; //clear any previous values - dispersion.clear(); + dispersion = rts::function< T, complex >(); //load the file into a string std::ifstream ifs(filename.c_str()); @@ -437,7 +441,6 @@ namespace rts{ //process the file as a string std::string instr((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); fromStr(instr, format, scaleA); - } void fromStr(std::string str, std::string format = "", T scaleA = 1.0) @@ -480,7 +483,6 @@ namespace rts{ //std::cout<<"Loading material with format: "<::compare); + //sort(dispersion.begin(), dispersion.end(), &material::compare); } //convert the material to a string @@ -555,51 +557,30 @@ namespace rts{ } - rtsComplex getN(T l) + complex getN(T l) { - //declare an iterator - typename std::vector< refIndex >::iterator it; - - refIndex r; - r.lambda = l; - - it = search(dispersion.begin(), dispersion.end(), &r, &r + 1, &material::findCeiling); - - //if the wavelength is past the end of the list, return the back - if(it == dispersion.end()) - return dispersion.back().n; - //if the wavelength is before the beginning of the list, return the front - else if(it == dispersion.begin()) - return dispersion.front().n; - //otherwise interpolate - else - { - T lMax = (*it).lambda; - T lMin = (*(it - 1)).lambda; - //std::cout< riMin = (*(it - 1)).n; - rtsComplex riMax = (*it).n; - rtsComplex interp; - interp = rtsComplex(a, 0.0) * riMin + rtsComplex(1.0 - a, 0.0) * riMax; - return interp; - } + //return complex(1.0, 0.0); + complex ri = dispersion.linear(l) + n0; + return ri; + } + function > getF() + { + return dispersion + complex(n0, 0.0); } //returns the scattering efficiency at wavelength l - rtsComplex eta(T l) + complex eta(T l) { //get the refractive index - rtsComplex ri = getN(l); + complex ri = getN(l); //convert the refractive index to scattering efficiency return ri*ri - 1.0; } //interpolate the given lambda value and return the index of refraction - rtsComplex operator()(T l) + complex operator()(T l) { return getN(l); } diff --git a/rts/source/colormap.cpp b/rts/source/colormap.cpp new file mode 100644 index 0000000..1176e26 --- /dev/null +++ b/rts/source/colormap.cpp @@ -0,0 +1,37 @@ +#include +#include +#include + +void qt_buffer2image(unsigned char* buffer, std::string filename, unsigned int x_size, unsigned int y_size) +{ + //x_size = 256; + //y_size = 256; + //create an image object + QImage image(x_size, y_size, QImage::Format_RGB32); + if(image.isNull()) + { + std::cout<<"Error creating QImage."< +#include +#include +#include + +#include + +#ifndef RTS_CAMERA_H +#define RTS_CAMERA_H + +namespace rts{ + +class camera +{ + vector d; //direction that the camera is pointing + point p; //position of the camera + vector up; //"up" direction + float focus; //focal length of the camera + float fov; + + //private function makes sure that the up vector is orthogonal to the direction vector and both are normalized + void stabalize() + { + vector side = up.cross(d); + up = d.cross(side); + up = up.norm(); + d = d.norm(); + } + +public: + void setPosition(point pos) + { + p = pos; + } + void setPosition(float x, float y, float z){setPosition(point(x, y, z));} + + void setFocalDistance(float distance){focus = distance;} + void setFOV(float field_of_view){fov = field_of_view;} + + void LookAt(point pos) + { + //find the new direction + d = pos - p; + + //find the distance from the look-at point to the current position + focus = d.len(); + + //stabalize the camera + stabalize(); + } + void LookAt(float px, float py, float pz){LookAt(point(px, py, pz));} + void LookAt(point pos, vector new_up){up = new_up; LookAt(pos);} + void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(point(px, py, pz), vector(ux, uy, uz));} + void LookAtDolly(float lx, float ly, float lz) + { + //find the current focus point + point f = p + focus*d; + vector T = point(lx, ly, lz) - f; + p = p + T; + } + + void Dolly(vector direction) + { + p = p+direction; + } + void Dolly(float x, float y, float z){Dolly(vector(x, y, z));} + void Push(float delta) + { + if(delta > focus) + delta = focus; + + focus -= delta; + + Dolly(d*delta); + } + + void Pan(float theta_x, float theta_y, float theta_z) + { + //x rotation is around the up axis + quaternion qx; + qx.CreateRotation(theta_x, up[0], up[1], up[2]); + + //y rotation is around the side axis + vector side = up.cross(d); + quaternion qy; + qy.CreateRotation(theta_y, side[0], side[1], side[2]); + + //z rotation is around the direction vector + quaternion qz; + qz.CreateRotation(theta_z, d[0], d[1], d[2]); + + //combine the rotations in x, y, z order + quaternion final = qz*qy*qx; + + //get the rotation matrix + matrix rot_matrix = final.toMatrix3(); + + //apply the rotation + d = rot_matrix*d; + up = rot_matrix*up; + + //stabalize the camera + stabalize(); + + } + void Pan(float theta_x){Pan(theta_x, 0, 0);} + void Tilt(float theta_y){Pan(0, theta_y, 0);} + void Twist(float theta_z){Pan(0, 0, theta_z);} + + void Zoom(float delta) + { + fov -= delta; + if(fov < 0.5) + fov = 0.5; + if(fov > 180) + fov = 180; + } + + void OrbitFocus(float theta_x, float theta_y) + { + //find the focal point + point focal_point = p + focus*d; + + //center the coordinate system on the focal point + point centered = p - (focal_point - point(0, 0, 0)); + + //create the x rotation (around the up vector) + quaternion qx; + qx.CreateRotation(theta_x, up[0], up[1], up[2]); + centered = point(0, 0, 0) + qx.toMatrix3()*(centered - point(0, 0, 0)); + + //get a side vector for theta_y rotation + vector side = up.cross((point(0, 0, 0) - centered).norm()); + + quaternion qy; + qy.CreateRotation(theta_y, side[0], side[1], side[2]); + centered = point(0, 0, 0) + qy.toMatrix3()*(centered - point(0, 0, 0)); + + //perform the rotation on the centered camera position + //centered = final.toMatrix()*centered; + + //re-position the camera + p = centered + (focal_point - point(0, 0, 0)); + + //make sure we are looking at the focal point + LookAt(focal_point); + + //stabalize the camera + stabalize(); + + } + + void Slide(float u, float v) + { + vector V = up.norm(); + vector U = up.cross(d).norm(); + + p = p + (V * v) + (U * u); + } + + //accessor methods + point getPosition(){return p;} + vector getUp(){return up;} + vector getDirection(){return d;} + point getLookAt(){return p + focus*d;} + float getFOV(){return fov;} + + //output the camera settings + const void print(std::ostream& output) + { + output<<"Position: "<(0, 0, 0); + d = vector(0, 0, 1); + up = vector(0, 1, 0); + focus = 1; + + } +}; + +} + + + +#endif diff --git a/rts/visualization/colormap.h b/rts/visualization/colormap.h new file mode 100644 index 0000000..c0be6f8 --- /dev/null +++ b/rts/visualization/colormap.h @@ -0,0 +1,321 @@ +#ifndef RTS_COLORMAP_H +#define RTS_COLORMAP_H + +#include +#include +#include "rts/cuda/error.h" + + +#define BREWER_CTRL_PTS 11 + +void qt_buffer2image(unsigned char* buffer, std::string filename, unsigned int x_size, unsigned int y_size); + +static float BREWERCP[BREWER_CTRL_PTS*4] = {0.192157f, 0.211765f, 0.584314f, 1.0f, + 0.270588f, 0.458824f, 0.705882f, 1.0f, + 0.454902f, 0.678431f, 0.819608f, 1.0f, + 0.670588f, 0.85098f, 0.913725f, 1.0f, + 0.878431f, 0.952941f, 0.972549f, 1.0f, + 1.0f, 1.0f, 0.74902f, 1.0f, + 0.996078f, 0.878431f, 0.564706f, 1.0f, + 0.992157f, 0.682353f, 0.380392f, 1.0f, + 0.956863f, 0.427451f, 0.262745f, 1.0f, + 0.843137f, 0.188235f, 0.152941f, 1.0f, + 0.647059f, 0.0f, 0.14902f, 1.0f}; + + +#ifdef __CUDACC__ +texture cudaTexBrewer; +static cudaArray* gpuBrewer; +#endif + + + +namespace rts{ + +enum colormapType {cmBrewer, cmGrayscale}; + +static void buffer2image(unsigned char* buffer, std::string filename, unsigned int x_size, unsigned int y_size) +{ + qt_buffer2image(buffer, filename, x_size, y_size); +} + +#ifdef __CUDACC__ +static void initBrewer() +{ + //initialize the Brewer colormap + + //allocate CPU space + float4 cpuColorMap[BREWER_CTRL_PTS]; + + //define control rtsPoints + cpuColorMap[0] = make_float4(0.192157f, 0.211765f, 0.584314f, 1.0f); + cpuColorMap[1] = make_float4(0.270588f, 0.458824f, 0.705882f, 1.0f); + cpuColorMap[2] = make_float4(0.454902f, 0.678431f, 0.819608f, 1.0f); + cpuColorMap[3] = make_float4(0.670588f, 0.85098f, 0.913725f, 1.0f); + cpuColorMap[4] = make_float4(0.878431f, 0.952941f, 0.972549f, 1.0f); + cpuColorMap[5] = make_float4(1.0f, 1.0f, 0.74902f, 1.0f); + cpuColorMap[6] = make_float4(0.996078f, 0.878431f, 0.564706f, 1.0f); + cpuColorMap[7] = make_float4(0.992157f, 0.682353f, 0.380392f, 1.0f); + cpuColorMap[8] = make_float4(0.956863f, 0.427451f, 0.262745f, 1.0f); + cpuColorMap[9] = make_float4(0.843137f, 0.188235f, 0.152941f, 1.0f); + cpuColorMap[10] = make_float4(0.647059f, 0.0f, 0.14902f, 1.0f); + + + int width = BREWER_CTRL_PTS; + int height = 0; + + + // allocate array and copy colormap data + cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat); + + HANDLE_ERROR(cudaMallocArray(&gpuBrewer, &channelDesc, width, height)); + + HANDLE_ERROR(cudaMemcpyToArray(gpuBrewer, 0, 0, cpuColorMap, sizeof(float4)*width, cudaMemcpyHostToDevice)); + + // set texture parameters + cudaTexBrewer.addressMode[0] = cudaAddressModeClamp; + //texBrewer.addressMode[1] = cudaAddressModeClamp; + cudaTexBrewer.filterMode = cudaFilterModeLinear; + cudaTexBrewer.normalized = true; // access with normalized texture coordinates + + // Bind the array to the texture + HANDLE_ERROR(cudaBindTextureToArray( cudaTexBrewer, gpuBrewer, channelDesc)); + +} + +static void destroyBrewer() +{ + HANDLE_ERROR(cudaFreeArray(gpuBrewer)); + +} + +template +__global__ static void applyBrewer(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1) +{ + + int i = blockIdx.y * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x; + if(i >= N) return; + + //compute the normalized value on [minVal maxVal] + float a = (gpuSource[i] - minVal) / (maxVal - minVal); + + //lookup the color + float shift = 1.0/(2*BREWER_CTRL_PTS); + float4 color = tex1D(cudaTexBrewer, a+shift); + //float4 color = tex1D(cudaTexBrewer, a); + + gpuDest[i * 3 + 0] = 255 * color.x; + gpuDest[i * 3 + 1] = 255 * color.y; + gpuDest[i * 3 + 2] = 255 * color.z; +} + +template +__global__ static void applyGrayscale(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1) +{ + int i = blockIdx.y * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x; + if(i >= N) return; + + //compute the normalized value on [minVal maxVal] + float a = (gpuSource[i] - minVal) / (maxVal - minVal); + + //threshold + if(a > 1.0) + a = 1.0; + if(a < 0.0) + a = 0.0; + + gpuDest[i * 3 + 0] = 255 * a; + gpuDest[i * 3 + 1] = 255 * a; + gpuDest[i * 3 + 2] = 255 * a; +} + +template +static void gpu2gpu(T* gpuSource, unsigned char* gpuDest, unsigned int nVals, T minVal = 0, T maxVal = 1, colormapType cm = cmGrayscale, int blockDim = 128) +{ + //This function converts a scalar field on the GPU to a color image on the GPU + int gridX = (nVals + blockDim - 1)/blockDim; + int gridY = 1; + if(gridX > 65535) + { + gridY = (gridX + 65535 - 1) / 65535; + gridX = 65535; + } + dim3 dimGrid(gridX, gridY); + //int gridDim = (nVals + blockDim - 1)/blockDim; + if(cm == cmGrayscale) + applyGrayscale<<>>(gpuSource, gpuDest, nVals, minVal, maxVal); + else if(cm == cmBrewer) + { + initBrewer(); + applyBrewer<<>>(gpuSource, gpuDest, nVals, minVal, maxVal); + //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3)); + destroyBrewer(); + } + +} + +template +static void gpu2cpu(T* gpuSource, unsigned char* cpuDest, unsigned int nVals, T minVal, T maxVal, colormapType cm = cmGrayscale) +{ + //this function converts a scalar field on the GPU to a color image on the CPU + + //first create the color image on the GPU + + //allocate GPU memory for the color image + unsigned char* gpuDest; + HANDLE_ERROR(cudaMalloc( (void**)&gpuDest, sizeof(unsigned char) * nVals * 3 )); + + //HANDLE_ERROR(cudaMemset(gpuSource, 0, sizeof(T) * nVals)); + + //create the image on the gpu + gpu2gpu(gpuSource, gpuDest, nVals, minVal, maxVal, cm); + + //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3)); + + //copy the image from the GPU to the CPU + HANDLE_ERROR(cudaMemcpy(cpuDest, gpuDest, sizeof(unsigned char) * nVals * 3, cudaMemcpyDeviceToHost)); + + HANDLE_ERROR(cudaFree( gpuDest )); + +} + +template +static void gpu2image(T* gpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale) +{ + //allocate a color buffer + unsigned char* cpuBuffer = NULL; + cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); + + //do the mapping + gpu2cpu(gpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm); + + //copy the buffer to an image + buffer2image(cpuBuffer, fileDest, x_size, y_size); + + free(cpuBuffer); +} + +#endif + +template +static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, unsigned int N, T minVal = 0, T maxVal = 1) +{ + for(int i=0; i 1) a = 1; + + float c = a * (float)(BREWER_CTRL_PTS-1); + int ptLow = (int)c; + float m = c - (float)ptLow; + //std::cout< +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T valMin, T valMax, colormapType cm = cmGrayscale) +{ + + if(cm == cmBrewer) + cpuApplyBrewer(cpuSource, cpuDest, nVals, valMin, valMax); + else if(cm == cmGrayscale) + { + int i; + float a; + float range = valMax - valMin; + for(i = 0; i 1) a = 1.0; + + cpuDest[i * 3 + 0] = 255 * a; + cpuDest[i * 3 + 1] = 255 * a; + cpuDest[i * 3 + 2] = 255 * a; + } + } +} + +template +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale, bool positive = false) +{ + //computes the max and min range automatically + + //find the largest magnitude value + T maxVal = fabs(cpuSource[0]); + for(int i=0; i maxVal) + maxVal = fabs(cpuSource[i]); + } + + if(positive) + cpu2cpu(cpuSource, cpuDest, nVals, (T)0.0, maxVal, cm); + else + cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm); + +} + + + +template +static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale) +{ + //allocate a color buffer + unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); + + //do the mapping + cpu2cpu(cpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm); + + //copy the buffer to an image + buffer2image(cpuBuffer, fileDest, x_size, y_size); + + free(cpuBuffer); + +} + +template +static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale, bool positive = false) +{ + //allocate a color buffer + unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size); + + //do the mapping + cpu2cpu(cpuSource, cpuBuffer, x_size * y_size, cm, positive); + + //copy the buffer to an image + buffer2image(cpuBuffer, fileDest, x_size, y_size); + + free(cpuBuffer); + +} + +} //end namespace colormap and rts + +#endif + -- libgit2 0.21.4