GPU单通道光线投射空体积（OpenGL，GLSL）

Question

你好，下面的代码是我的一个简单体积单通道光线投射算法的渲染管道。它使用 OpenGL (glew,freeglut) 和 GLSL。当我渲染体积时，我得到的只是体积本身的外部，但如果我滚动到它内部则不是。我会提供截图。 Shaderloader 非常简单，所以我不会在这里发布。这是我的代码（我删除了一些回调以使问题更短）：

//OpenGL
#include <GL/glew.h>
#include <GL/freeglut.h>

//Base C++
#include <fstream>
#include <iostream>
#include <chrono>

//GLM
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>

//Shaderloader
#include "GLSLShader.h"


#define GL_CHECK_ERRORS assert(glGetError()== GL_NO_ERROR);


using namespace std;

//main screen resolution
const int WIDTH = 1280;
const int HEIGHT = 960;

//second screen resolution 
const int WIDTH2 = 640;
const int HEIGHT2 = 480;
const int POS_X2 = 50;
const int POS_Y2 = 50;

//camera transform variables
int state = 0, oldX = 0, oldY = 0;
float rX = 4, rY = 50, dist = -2;

//modelview projection matrices
glm::mat4 MV, P;

//cube vertex array and vertex buffer object IDs
GLuint cubeVBOID;
GLuint cubeVAOID;
GLuint cubeIndicesID;

//ray casting shader
GLSLShader shader;

//background colour
glm::vec4 bg = glm::vec4(0.2, 0.2, 0.2, 1);

//volume dataset filename  
const std::string volume_file = "Engine256.raw";

//volume dimensions
const int XDIM = 256;
const int YDIM = 256;
const int ZDIM = 256;

//fps tracker
static int fps = 0;
static int frameCount = 0;
static int currentTime = 0;
static int previousTime = 0;

//volume texture ID
GLuint textureID;

//transfer function (lookup table) texture id
GLuint tfTexID;



//transfer function (lookup table) colour values
const glm::vec4 jet_values[9] = { glm::vec4(0,0,0.5,0),
                                glm::vec4(0,0,1,0.1),
                                glm::vec4(0,0.5,1,0.3),
                                glm::vec4(0,1,1,0.5),
                                glm::vec4(0.5,1,0.5,0.75),
                                glm::vec4(1,1,0,0.8),
                                glm::vec4(1,0.5,0,0.6),
                                glm::vec4(1,0,0,0.5),
                                glm::vec4(0.5,0,0,0.0) };

//function that load a volume from the given raw data file and 
//generates an OpenGL 3D texture from it
bool LoadVolume() {
    std::ifstream infile(volume_file.c_str(), std::ios_base::binary);

    if (infile.good()) {
        //read the volume data file
        GLubyte* pData = new GLubyte[XDIM * YDIM * ZDIM];
        infile.read(reinterpret_cast<char*>(pData), XDIM * YDIM * ZDIM * sizeof(GLubyte));
        infile.close();

        //generate OpenGL texture
        glGenTextures(1, &textureID);
        glBindTexture(GL_TEXTURE_3D, textureID);

        // set the texture parameters
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP);
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP);
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP);
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR);

        //set the mipmap levels (base and max)
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_BASE_LEVEL, 0);
        glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAX_LEVEL, 4);

        //allocate data with internal format and foramt as (GL_RED) 
        glTexImage3D(GL_TEXTURE_3D, 0, GL_RED, XDIM, YDIM, ZDIM, 0, GL_RED, GL_UNSIGNED_BYTE, pData);
        GL_CHECK_ERRORS

        //generate mipmaps
        glGenerateMipmap(GL_TEXTURE_3D);

        //delete the volume data allocated on heap
        delete[] pData;

        return true;
    }
    else {
        return false;
    }
}

//function to generate interpolated colours from the set of colour values (jet_values)
//this function first calculates the amount of increments for each component and the
//index difference. Then it linearly interpolates the adjacent values to get the 
//interpolated result.
void LoadTransferFunction() {
    float pData[256][4];
    int indices[9];

    //fill the colour values at the place where the colour should be after interpolation
    for (int i = 0; i < 9; i++) {
        int index = i * 28;
        pData[index][0] = jet_values[i].x;
        pData[index][1] = jet_values[i].y;
        pData[index][2] = jet_values[i].z;
        pData[index][3] = jet_values[i].w;
        indices[i] = index;
    }

    //for each adjacent pair of colours, find the difference in the rgba values and then interpolate
    for (int j = 0; j < 9 - 1; j++)
    {
        float dDataR = (pData[indices[j + 1]][0] - pData[indices[j]][0]);
        float dDataG = (pData[indices[j + 1]][1] - pData[indices[j]][1]);
        float dDataB = (pData[indices[j + 1]][2] - pData[indices[j]][2]);
        float dDataA = (pData[indices[j + 1]][3] - pData[indices[j]][3]);
        int dIndex = indices[j + 1] - indices[j];

        float dDataIncR = dDataR / float(dIndex);
        float dDataIncG = dDataG / float(dIndex);
        float dDataIncB = dDataB / float(dIndex);
        float dDataIncA = dDataA / float(dIndex);
        for (int i = indices[j] + 1; i < indices[j + 1]; i++)
        {
            pData[i][0] = (pData[i - 1][0] + dDataIncR);
            pData[i][1] = (pData[i - 1][1] + dDataIncG);
            pData[i][2] = (pData[i - 1][2] + dDataIncB);
            pData[i][3] = (pData[i - 1][3] + dDataIncA);
        }
    }

    //generate the OpenGL texture
    glGenTextures(1, &tfTexID);
    //bind this texture to texture unit 1
    glActiveTexture(GL_TEXTURE1);
    glBindTexture(GL_TEXTURE_1D, tfTexID);

    // set the texture parameters
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_REPEAT);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);

    //allocate the data to texture memory. Since pData is on stack, we donot delete it 
    glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA, 256, 0, GL_RGBA, GL_FLOAT, pData);

    GL_CHECK_ERRORS
}


//mouse down event handler
void OnMouseDown(int button, int s, int x, int y)

//mouse move event handler
void OnMouseMove(int x, int y)

//scroll callback function
void OnMouseScroll(int wheel, int direction, int x, int y)

//OpenGL initialization
void OnInit() {

    auto start = std::chrono::steady_clock::now();

    GL_CHECK_ERRORS

    //Load the raycasting shader
    shader.LoadFromFile(GL_VERTEX_SHADER, "raycaster.vert");
    shader.LoadFromFile(GL_FRAGMENT_SHADER, "raycaster.frag");


    //compile and link the shader
    shader.CreateAndLinkProgram();
    shader.Use();
    //add attributes and uniforms
    shader.AddAttribute("vVertex");
    shader.AddUniform("MVP");
    shader.AddUniform("volume");
    shader.AddUniform("camPos");
    shader.AddUniform("step_size");
    shader.AddUniform("lut");
    

    //pass constant uniforms at initialization
    glUniform3f(shader("step_size"), 1.0f / XDIM, 1.0f / YDIM, 1.0f / ZDIM);
    glUniform1i(shader("volume"), 0);
    glUniform1i(shader("lut"), 1);
    shader.UnUse();

    //GL_CHECK_ERRORS

        //load volume data
        if (LoadVolume()) {
            std::cout << "Volume data loaded successfully." << std::endl;
        }
        else {
            std::cout << "Cannot load volume data." << std::endl;
            exit(EXIT_FAILURE);
        }

//load the transfer function data and generate the trasnfer function (lookup table) texture
    LoadTransferFunction();

    //set background colour
    glClearColor(bg.r, bg.g, bg.b, bg.a);

    //setup unit cube vertex array and vertex buffer objects
    glGenVertexArrays(1, &cubeVAOID);
    glGenBuffers(1, &cubeVBOID);
    glGenBuffers(1, &cubeIndicesID);

    //unit cube vertices 
    glm::vec3 vertices[8] = { glm::vec3(-0.5f,-0.5f,-0.5f),
                            glm::vec3(0.5f,-0.5f,-0.5f),
                            glm::vec3(0.5f, 0.5f,-0.5f),
                            glm::vec3(-0.5f, 0.5f,-0.5f),
                            glm::vec3(-0.5f,-0.5f, 0.5f),
                            glm::vec3(0.5f,-0.5f, 0.5f),
                            glm::vec3(0.5f, 0.5f, 0.5f),
                            glm::vec3(-0.5f, 0.5f, 0.5f) };

    //unit cube indices
    GLushort cubeIndices[36] = { 0,5,4,
                              5,0,1,
                              3,7,6,
                              3,6,2,
                              7,4,6,
                              6,4,5,
                              2,1,3,
                              3,1,0,
                              3,0,7,
                              7,0,4,
                              6,5,2,
                              2,5,1 };
    glBindVertexArray(cubeVAOID);
    glBindBuffer(GL_ARRAY_BUFFER, cubeVBOID);
    //pass cube vertices to buffer object memory
    glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), &(vertices[0].x), GL_STATIC_DRAW);

    GL_CHECK_ERRORS
    //enable vertex attributre array for position
    glEnableVertexAttribArray(0);
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, 0);

    //pass indices to element array  buffer
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, cubeIndicesID);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(cubeIndices), &cubeIndices[0], GL_STATIC_DRAW);

    glBindVertexArray(0);

    //enable depth test
    glEnable(GL_DEPTH_TEST);

    //set the over blending function
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    cout << "Initialization successfull" << endl;
    auto end = std::chrono::steady_clock::now();
    chrono::duration<double> elapsed_seconds = end - start;
    cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
}

//release all allocated resources
void OnShutdown() {

//resize event handler
void OnResize(int w, int h) {

//display callback function
void OnRender() {
    GL_CHECK_ERRORS

    //set the camera transform
    glm::mat4 Tr = glm::translate(glm::mat4(1.0f), glm::vec3(0.0f, 0.0f, dist));
    glm::mat4 Rx = glm::rotate(Tr, rX, glm::vec3(1.0f, 0.0f, 0.0f));
    glm::mat4 MV = glm::rotate(Rx, rY, glm::vec3(0.0f, 1.0f, 0.0f));

    //get the camera position
    glm::vec3 camPos = glm::vec3(glm::inverse(MV) * glm::vec4(0, 0, 0, 1));

    //clear colour and depth buffer
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    //get the combined modelview projection matrix
    glm::mat4 MVP = P * MV;

    //enable blending and bind the cube vertex array object
    glEnable(GL_BLEND);
    glBindVertexArray(cubeVAOID);
    //bind the raycasting shader
    shader.Use();
    //pass shader uniforms
    glUniformMatrix4fv(shader("MVP"), 1, GL_FALSE, glm::value_ptr(MVP));
    glUniform3fv(shader("camPos"), 1, &(camPos.x));
    //render the cube
    glDrawElements(GL_TRIANGLES, 36, GL_UNSIGNED_SHORT, 0);
    //unbind the raycasting shader
    shader.UnUse();
    //disable blending
    glDisable(GL_BLEND);

    //swap front and back buffers to show the rendered result
    glutSwapBuffers();
}

int main(int argc, char** argv) {
    //freeglut initialization
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
    glutInitContextVersion(3, 3);
    glutInitContextFlags(GLUT_CORE_PROFILE | GLUT_DEBUG);

    //main window
    glutInitWindowSize(WIDTH, HEIGHT);
    glutCreateWindow("Volume Rendering using GPU Ray Casting - OpenGL 3.3");


    //glew initialization
    glewExperimental = GL_TRUE;
    GLenum err = glewInit();
    if (GLEW_OK != err) {
        cerr << "Error: " << glewGetErrorString(err) << endl;
    }
    else {
        if (GLEW_VERSION_3_3)
        {
            cout << "Driver supports OpenGL 3.3\nDetails:" << endl;
        }
    }
    err = glGetError(); //this is to ignore INVALID ENUM error 1282
    GL_CHECK_ERRORS

    //output hardware information
    cout << "\tUsing GLEW " << glewGetString(GLEW_VERSION) << endl;
    cout << "\tVendor: " << glGetString(GL_VENDOR) << endl;
    cout << "\tRenderer: " << glGetString(GL_RENDERER) << endl;
    cout << "\tVersion: " << glGetString(GL_VERSION) << endl;
    cout << "\tGLSL: " << glGetString(GL_SHADING_LANGUAGE_VERSION) << endl;

    GL_CHECK_ERRORS

    //OpenGL initialization
    OnInit();

    //callback hooks
    glutCloseFunc(OnShutdown);
    glutDisplayFunc(OnRender);
    glutReshapeFunc(OnResize);
    glutMouseFunc(OnMouseDown);
    glutMotionFunc(OnMouseMove);
    glutMouseWheelFunc(OnMouseScroll);

    //main loop call
    glutMainLoop();

    return 0;
}

顶点着色器：

#version 330 core
  
layout(location = 0) in vec3 vVertex; //object space vertex position

//uniform
uniform mat4 MVP;   //combined modelview projection matrix

smooth out vec3 vUV; //3D texture coordinates for texture lookup in the fragment shader

void main()
{  
//get the clipspace position 
gl_Position = MVP*vec4(vVertex.xyz,1);

//get the 3D texture coordinates by adding (0.5,0.5,0.5) to the object space 
//vertex position. Since the unit cube is at origin (min: (-0.5,-0.5,-0.5) and max: (0.5,0.5,0.5))
//adding (0.5,0.5,0.5) to the unit cube object space position gives us values from (0,0,0) to 
//(1,1,1)
    vUV = vVertex + vec3(0.5);
}

片段着色器：

version 330 core

layout(location = 0) out vec4 vFragColor;   //fragment shader output

smooth in vec3 vUV;             //3D texture coordinates form vertex shader 
                        //interpolated by rasterizer

//uniforms
uniform sampler3D   volume;     //volume dataset
uniform sampler1D   lut;        //look-up-table for transferfunction
uniform vec3        camPos;     //camera position
uniform vec3        step_size;  //ray step size 

//constants
const int MAX_SAMPLES = 300;    //total samples for each ray march step
const vec3 texMin = vec3(0);    //minimum texture access coordinate
const vec3 texMax = vec3(1);    //maximum texture access coordinate

void main()
{ 
    //transferfunction
    vFragColor = texture(lut, texture(volume, vUV).r);

    //get the 3D texture coordinates for lookup into the volume dataset
    vec3 dataPos = vUV;

    //Getting the ray marching direction:
    //get the object space position by subracting 0.5 from the
    //3D texture coordinates. Then subtraact it from camera position
    //and normalize to get the ray marching direction
    vec3 geomDir = normalize((vUV-vec3(0.5)) - camPos); 

    //multiply the raymarching direction with the step size to get the
    //sub-step size we need to take at each raymarching step
    vec3 dirStep = geomDir * step_size; 
     
    //flag to indicate if the raymarch loop should terminate
    bool stop = false; 

    //for all samples along the ray
    for (int i = 0; i < MAX_SAMPLES; i++) {
        // advance ray by dirstep
        dataPos = dataPos + dirStep;
        
        
        //The two constants texMin and texMax have a value of vec3(-1,-1,-1)
        //and vec3(1,1,1) respectively. To determine if the data value is 
        //outside the volume data, we use the sign function. The sign function 
        //return -1 if the value is less than 0, 0 if the value is equal to 0 
        //and 1 if value is greater than 0. Hence, the sign function for the 
        //calculation (sign(dataPos-texMin) and sign (texMax-dataPos)) will 
        //give us vec3(1,1,1) at the possible minimum and maximum position. 
        //When we do a dot product between two vec3(1,1,1) we get the answer 3. 
        //So to be within the dataset limits, the dot product will return a 
        //value less than 3. If it is greater than 3, we are already out of 
        //the volume dataset
        stop = dot(sign(dataPos-texMin),sign(texMax-dataPos)) < 3.0;

        //if the stopping condition is true we brek out of the ray marching loop
        if (stop) 
            break;

        // data fetching from the red channel of volume texture
        float sample = texture(volume, dataPos).r;  
        
        //Opacity calculation using compositing:
        //here we use front to back compositing scheme whereby the current sample
        //value is multiplied to the currently accumulated alpha and then this product
        //is subtracted from the sample value to get the alpha from the previous steps.
        //Next, this alpha is multiplied with the current sample colour and accumulated
        //to the composited colour. The alpha value from the previous steps is then 
        //accumulated to the composited colour alpha.
        float prev_alpha = sample - (sample * vFragColor.a);
        vFragColor.rgb = prev_alpha * vec3(sample) + vFragColor.rgb; 
        vFragColor.a += prev_alpha; 
            
        //early ray termination
        //if the currently composited colour alpha is already fully saturated
        //we terminated the loop
        if( vFragColor.a>0.99)
            break;
    } 
}