gpumathlib.cpp

这是在GPU上进行高速数字计算的数学库函数

/*******************************************************************************
**
** The header file for the gpumathlib class
** The Jahshaka Project
** Copyright (C) 2000-2004 The Jahshaka Project.
** Released under the GNU General Public License
**
*******************************************************************************/

#include "gpumathlib.h"
#include <assert.h>

// Utility routines here
float sfrand(void)
{
    return (rand() * 2.0f / (float)RAND_MAX) - 1.0f;
}

float frand(void)
{
    return (rand() / (float)RAND_MAX);
}

void normalize3f(float& x, float& y, float& z)
{
    float						length;

    length = sqrt((x * x) + (y * y) + (z * z));

    x /= length;
    y /= length;
    z /= length;
}

float normalizefloat4(float4& vector)
{
    float						length;

    length = sqrt((vector.x * vector.x) + (vector.y * vector.y) + (vector.z * vector.z) + (vector.w * vector.w));

    if (length > .000001f)
    {
        vector.x /= length;
        vector.y /= length;
        vector.z /= length;
        vector.w /= length;
    }

    return length;
}

void get_cross_product_4( float4& result, float4 vector1, float4 vector2)
{
    // Assumes that w is 1
    result.x = (vector1.y * vector2.z) - (vector1.z * vector2.y);
    result.y = -((vector1.x * vector2.z) - (vector1.z * vector2.x));
    result.z = (vector1.x * vector2.y) - (vector1.y * vector2.x);
    result.w = 1.0;
}

void normalize_vector4(float4& vector)
{
    float					length;

    length = sqrt( (vector.x * vector.x) + (vector.y * vector.y) + (vector.z * vector.z) );

    vector.x /= length;
    vector.y /= length;
    vector.z /= length;
    vector.w = 1.0;
// 	vector.w /= length;
}

void get_vector_from_points(float4& vector, float4 point1, float4 point2)
{

    vector.x = point2.x - point1.x;
    vector.y = point2.y - point1.y;
    vector.z = point2.z - point1.z;
    vector.w = 1.0;
// 	vector.w = point2.w - point1.w;

    normalize_vector4(vector);
}


static char error_txt[][32] =
{
    "GL_INVALID_ENUM",
    "GL_INVALID_VALUE",
    "GL_INVALID_OPERATION",
    "GL_STACK_OVERFLOW",
    "GL_STACK_UNDERFLOW",
    "GL_OUT_OF_MEMORY"
};

static const unsigned int max_errors = 6;


void check_gl()
{
    GLenum                gl_error;
    int                    error_offset;

    gl_error = glGetError();
    error_offset = gl_error - GL_INVALID_ENUM;


    if (gl_error != GL_NO_ERROR)
    {
        if (error_offset >= (int)max_errors)
        {
            fprintf(stderr, "Unknown GL error %08x\n", gl_error);
        }
        else
        {
            fprintf(stderr, "GL ERROR %s\n", error_txt[error_offset]);
        }

        assert(0);
    }
}

///////////////////////////////////////////////////
// this crashes when returning NULL
// we need to pass the return value over after name
// and return a bool if it worked or not

unsigned char* loadshaderfile(const char *name)
{
    FILE*                       file;
    int                         size;
    unsigned char*              data;

    file = fopen(name, "rb");

    if(!file)
    {
        fprintf(stderr, "loadshaderfile:: error loading %s file",name);
        return NULL; 
    }

    fseek(file,0,SEEK_END);
    size = ftell(file);
    data = new unsigned char[size + 1];
    fseek(file,0,SEEK_SET);
    fread(data,1,size,file);
    data[size] = '\0';
    fclose(file);
    return data;
}



void find_shader_program_error(unsigned char * source, char* program_name)
{

    GLint                           pos;
    int                             i;
    int                             line;
    int                             linestart;
    char*                           progcopy;


    progcopy = strdup((char*)source);
    glGetIntegerv(GL_PROGRAM_ERROR_POSITION_NV, &pos);

    line = 1;
    linestart = 0;

    for (i = 0; i < pos; i++)
    {
        if (progcopy[i] == '\n')
        {
            line++;
            linestart = i+1;
        }
    }

    if (line == 1)
    {
        return;
    }

    fprintf ( stderr, "find_shader_program_error: Program Error : %s : line %d\n",
            program_name, line);

    for (i = linestart; progcopy[i] != '\0' && progcopy[i] != '\n'; i++);

    progcopy[i] = '\0';

    fprintf ( stderr, "%s\n", progcopy+linestart);

    for (i=linestart; i<pos; i++)
    {
        fprintf ( stderr, " ");
    }

    for (;progcopy[i] != '\0' && progcopy[i] != '\n'; i++) 
    {
        fprintf ( stderr, "^");
    }

    fprintf ( stderr, "\n");
    free(progcopy);
    fprintf ( stderr, "%s\n", glGetString(GL_PROGRAM_ERROR_STRING_ARB));

}

void 
matrixTranspose4x4(float* matrix)
{
    float temp[4][4];
    int row;
    int column;

    for (row = 0; row < 4; row++)
    {
        for (column = 0; column < 4; column++)
        {
            temp[row][column] = matrix[(column * 4) + row];
        }
    }

    for (row = 0; row < 4; row++)
    {
        for (column = 0; column < 4; column++)
        {
            matrix[(row * 4) + column] = temp[row][column];
        }
    }
}

void
matrixMultiply4x4(float* matrix_out, float* matrix1, float* matrix2)
{
    int         row;
    int         column;

    for (row = 0; row < 4; row++)
    {
        for (column = 0; column < 4; column++)
        {
            matrix_out[(row * 4) + column] = 0.0f;

            for (int index = 0; index < 4; index++)
            {
                matrix_out[(row * 4) + column] += 
                        matrix1[(row * 4) + index] * 
                        matrix2[(index * 4) + column];
            }
        }
    }
}

void 
getMVPMatrices(float* modelview, float* projection, float* mvp, float* modelviewIT, float* texture)
{
    glGetFloatv(GL_MODELVIEW_MATRIX, modelview);
    glGetFloatv(GL_PROJECTION_MATRIX, projection);
    glGetFloatv(GL_TEXTURE_MATRIX, texture);

    matrixMultiply4x4(mvp, modelview, projection);

    // Put them in row major to pass to the shader
    matrixTranspose4x4(modelview);
    matrixTranspose4x4(projection);
    matrixTranspose4x4(mvp);
    matrixTranspose4x4(texture);

    //matrixPrint4x4(modelview, "MODELVIEW");
    //matrixPrint4x4(projection, "PROJECTION");
    //matrixPrint4x4(texture, "texture");

    matrixInvert4x4(modelviewIT, modelview);
    // Already transposed by OpenGL?
    // matrixTranspose4x4(modelviewIT);
}

void
matrixPrintTranspose4x4(float* matrix, const char* name)
{
    fprintf(stderr, "%s :\n", name);

    for (int row = 0; row < 4; row++)
    {
        fprintf(stderr, "%8.3f, %8.3f, %8.3f, %8.3f\n", 
                matrix[0 + row],
                matrix[4 + row],
                matrix[8 + row],
                matrix[12 + row]);
    }
}

void
matrixPrint4x4(float* matrix, const char* name)
{
    fprintf(stderr, "%s :\n", name);

    for (int row = 0; row < 4; row++)
    {
        fprintf(stderr, "%8.3f, %8.3f, %8.3f, %8.3f\n", 
                matrix[(row * 4) + 0],
                matrix[(row * 4) + 1],
                matrix[(row * 4) + 2],
                matrix[(row * 4) + 3]);
    }
}

void
matrixSwapRows4x4(float* row1, float* row2)
{
    int column;
    float temp;

    for (column = 0; column < 4; column++)
    {
        temp = row1[column];
        row1[column] = row2[column];
        row2[column] = temp;
    }
}

void 
matrixInvert4x4(float* inverse, float* matrix)
{
    int         row;
    int         column;
    float       factor;
    int         diagonal;


    static const bool debug_level = false;

    // Using gauss-jordan elimination

    // Load the identity matrix
    for (row = 0; row < 4; row++)
    {
        for (column = 0; column < 4; column++)
        {
            if (row == column)
            {
                inverse[(row * 4) + column] = 1.0f;
            }
            else
            {
                inverse[(row * 4) + column] = 0.0f;
            }
        }
    }

    if (debug_level)
    {
        matrixPrint4x4(matrix, "Initial MATRIX");
        matrixPrint4x4(inverse, "Initial INVERSE");
    }

    for (diagonal = 0; diagonal < 4; diagonal++)
    {
        // If the diagonal element is 0, try to swap with another row from below
        if (matrix[(diagonal * 4) + diagonal] == 0.0f)
        {
            for (row = diagonal; row < 4; row++)
            {
                if (matrix[(row * 4) + diagonal] != 0.0f)
                {
                    matrixSwapRows4x4(matrix + (diagonal * 4), matrix + (row * 4));
                    continue;
                }

                debug("matrixInvert4x4:: Matrix is non-invertible\n");
                return;
            }
        }

        for (row = 0; row < 4; row++)
        {
            if (row != diagonal)
            {
                factor = (matrix[(row * 4) + diagonal] / matrix[(diagonal * 4) + diagonal]);
                
                for (column = 0; column < 4; column++)
                {
                    matrix[(row * 4) + column] -= (matrix[(diagonal * 4) + column] * factor);
                    inverse[(row * 4) + column] -= (inverse[(diagonal * 4) + column] * factor);
                }

                if (debug_level)
                {
                    matrixPrint4x4(matrix, "MATRIX row conversion");
                    matrixPrint4x4(inverse, "INVERSE row conversion");
                }