Volume Rendering

[Joaopmoliveira]

Good afternoon VSG community,

I have been trying to port the OpenSceneGraph volume rendering shaders into the VSG API, and something is sort of working.

I would like to get some help from people with more experience with the API to add this feature to VSG (if no help is possible, general directions of how to implement the feature are welcome).

#include <vsg/all.h>
#include <vsgXchange/all.h>
#include <iostream>

char volume_vert[] = R"(
#version 450
#extension GL_ARB_separate_shader_objects : enable

layout(push_constant) uniform PushConstants {
    mat4 projection;
    mat4 modelview;
} pc;

layout(location = 0) in vec3 inPosition;

layout(location = 0) out vec4 cameraPos;
layout(location = 1) out vec4 vertexPos;
layout(location = 2) out mat4 texgen;

out gl_PerVertex{
    vec4 gl_Position;
};

void main() {
    gl_Position = (pc.projection * pc.modelview) * vec4(inPosition, 1.0);
    cameraPos = inverse(pc.modelview) * vec4(0,0,0,1);
    vertexPos = vec4(inPosition, 1.0);
    texgen = mat4(1.0);
})";


char volume_frag[] = R"(
#version 450
#extension GL_ARB_separate_shader_objects : enable

layout(binding = 0) uniform sampler3D volume;

layout(location = 0) in vec4 cameraPos;
layout(location = 1) in vec4 vertexPos;
layout(location = 2) in mat4 texgen;

layout(location = 0) out vec4 outColor;

void main() {
    vec4 t0 = vertexPos;
    vec4 te = cameraPos;
    if( te.x>=0.0 && te.x<=1.0 &&
        te.y>=0.0 && te.y<=1.0 &&
        te.z>=0.0 && te.z<=1.0
        )
    { } else {
        if (te.x<0.0){
            float r = -te.x / (t0.x-te.x);
            te = te + (t0-te)*r;
        }
        if (te.x>1.0){
            float r = (1.0-te.x) / (t0.x-te.x);
            te = te + (t0-te)*r;
        }
        if (te.y<0.0){
            float r = -te.y / (t0.y-te.y);
            te = te + (t0-te)*r;
        }
        if (te.y>1.0){
            float r = (1.0-te.y) / (t0.y-te.y);
            te = te + (t0-te)*r;
        }
        if (te.z<0.0){
            float r = -te.z / (t0.z-te.z);
            te = te + (t0-te)*r;
        }
        if (te.z>1.0){
            float r = (1.0-te.z) / (t0.z-te.z);
            te = te + (t0-te)*r;
        }
    }

    t0 = t0 * texgen;
    te = te * texgen;

    const float min_iteratrions = 2.0;
    const float max_iteratrions = 2048.0;

    float TransparencyValue = 0.8;
    float AlphaFuncValue = 0.2;
    float SampleDensityValue = 10;
    
    float num_iterations = ceil(length((te-t0).xyz)/SampleDensityValue);
    if (num_iterations<min_iteratrions) num_iterations = min_iteratrions;
    else if (num_iterations>max_iteratrions) num_iterations = max_iteratrions;

    vec3 deltaTexCoord=(te-t0).xyz/(num_iterations-1.0);
    vec3 texcoord = t0.xyz;

    vec4 fragColor = vec4(0.0, 0.0, 0.0, 0.0);
    while(num_iterations>0.0){
        vec4 color1 = texture( volume, texcoord);
        vec4 color = vec4(1-color1.r,1-color1.r,1-color1.r,color1.r);
        float r = color[3]*TransparencyValue;
        if (r>AlphaFuncValue)
        {
            fragColor.xyz = fragColor.xyz*(1.0-r)+color.xyz*r;
            fragColor.w += r;
        }
        if (fragColor.w<color.w){
            fragColor = color;
        }
        texcoord += deltaTexCoord;
        --num_iterations;
    }
    fragColor.w *= TransparencyValue;
    if (fragColor.w>1.0) fragColor.w = 1.0;
    //if (fragColor.w<AlphaFuncValue) discard;
    outColor = fragColor;
})";

void updateBaseTexture3D(vsg::floatArray3D& image, float value)
{
    for (size_t d = 0; d < image.depth(); ++d) {
        float d_ratio = static_cast<float>(d) / static_cast<float>(image.depth() - 1);
        for (size_t r = 0; r < image.height(); ++r)
        {
            float r_ratio = static_cast<float>(r) / static_cast<float>(image.height() - 1);
            for (size_t c = 0; c < image.width(); ++c)
            {
                float c_ratio = static_cast<float>(c) / static_cast<float>(image.width() - 1);
                
                vsg::vec3 delta((r_ratio - 0.5f), (c_ratio - 0.5f), (d_ratio - 0.5));
                float angle = atan2(delta.x, delta.y);
                float distance_from_center = sqrt(d_ratio * d_ratio + r_ratio * r_ratio);

                float intensity = (sin(1.0 * angle + 30.0f * distance_from_center + 10.0 * value) + 1.0f) * 0.5f;
                image.set(c, r, d, 1.0f-intensity);
            }
        }
    }

    image.dirty();
}


vsg::ref_ptr<vsg::Node> get_volume_rendering() {
    vsg::vec3 origin = vsg::vec3(0.0, 0.0, 0.0);
    vsg::vec3 dx = vsg::vec3(1.0, 0.0, 0.0);
    vsg::vec3 dy = vsg::vec3(0.0, 1.0, 0.0);
    vsg::vec3 dz = vsg::vec3(0.0, 0.0, 1.0);
    vsg::vec3 v000(origin);
    vsg::vec3 v100(origin + dx);
    vsg::vec3 v110(origin + dx + dy);
    vsg::vec3 v010(origin + dy);
    vsg::vec3 v001(origin + dz);
    vsg::vec3 v101(origin + dx + dz);
    vsg::vec3 v111(origin + dx + dy + dz);
    vsg::vec3 v011(origin + dy + dz);

    vsg::ref_ptr<vsg::vec3Array> vertices;
    vsg::ref_ptr<vsg::vec3Array> normals;
    vsg::ref_ptr<vsg::vec2Array> texcoords;
    vsg::ref_ptr<vsg::ushortArray> indices;

    vsg::vec3 n0 = normalize(cross(dx, dz));
    vsg::vec3 n1 = normalize(cross(dy, dz));
    vsg::vec3 n2 = -n0;
    vsg::vec3 n3 = -n1;
    vsg::vec3 n4 = normalize(cross(dy, dx));
    vsg::vec3 n5 = -n4;

    // set up vertex and index arrays
    vertices = vsg::vec3Array::create(
        { v000, v100, v101, v001,   // front
         v100, v110, v111, v101,   // right
         v110, v010, v011, v111,   // far
         v010, v000, v001, v011,   // left
         v010, v110, v100, v000,   // bottom
         v001, v101, v111, v011 }); // top

    indices = vsg::ushortArray::create(
        { 0, 1, 2, 0, 2, 3,
         4, 5, 6, 4, 6, 7,
         8, 9, 10, 8, 10, 11,
         12, 13, 14, 12, 14, 15,
         16, 17, 18, 16, 18, 19,
         20, 21, 22, 20, 22, 23 });

    vsg::DataList arrays;
    arrays.push_back(vertices);
    arrays.push_back(normals);

    auto vid = vsg::VertexIndexDraw::create();
    vid->assignArrays(arrays);
    vid->assignIndices(indices);
    vid->indexCount = static_cast<uint32_t>(indices->size());
    uint32_t instanceCount = 1;
    vid->instanceCount = instanceCount;
    return vid;
};


int main(int argc, char** argv) {
    // set up defaults and read command line arguments to override them
    vsg::CommandLine arguments(&argc, argv);

    auto windowTraits = vsg::WindowTraits::create();
    windowTraits->debugLayer = true;
    windowTraits->apiDumpLayer = false;
    windowTraits->width = 1000;
    windowTraits->height = 1000;

    if (arguments.errors()) return arguments.writeErrorMessages(std::cerr);


    // load shaders
    auto vertexShader = vsg::ShaderStage::create(VK_SHADER_STAGE_VERTEX_BIT, "main", volume_vert);
    auto fragmentShader = vsg::ShaderStage::create(VK_SHADER_STAGE_FRAGMENT_BIT, "main", volume_frag);

    if (!vertexShader || !fragmentShader)
    {
        std::cout << "Could not create shaders." << std::endl;
        return {};
    }

    // read texture image
    auto textureData = vsg::floatArray3D::create(100, 100, 100);
    textureData->properties.format = VK_FORMAT_R32_SFLOAT;
    textureData->properties.dataVariance = vsg::DYNAMIC_DATA;

    updateBaseTexture3D(*textureData, 1.0f);
    if (!textureData)
    {
        std::cout << "Could not read texture" << std::endl;
        return {};
    }

    // set up graphics pipeline
    vsg::DescriptorSetLayoutBindings descriptorBindings{
        {0, VK_DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER, 1, VK_SHADER_STAGE_FRAGMENT_BIT, nullptr} // { binding, descriptorTpe, descriptorCount, stageFlags, pImmutableSamplers}
    };

    auto descriptorSetLayout = vsg::DescriptorSetLayout::create(descriptorBindings);

    vsg::PushConstantRanges pushConstantRanges{
        {VK_SHADER_STAGE_VERTEX_BIT, 0, 128} // projection view, and model matrices, actual push constant calls automatically provided by the VSG's DispatchTraversal
    };

    vsg::VertexInputState::Bindings vertexBindingsDescriptions{
        VkVertexInputBindingDescription{0, sizeof(vsg::vec3), VK_VERTEX_INPUT_RATE_VERTEX}, // vertex data
    };

    vsg::VertexInputState::Attributes vertexAttributeDescriptions{
        VkVertexInputAttributeDescription{0, 0, VK_FORMAT_R32G32B32_SFLOAT, 0}, // vertex data
    };

    vsg::GraphicsPipelineStates pipelineStates{
      vsg::VertexInputState::create(vertexBindingsDescriptions, vertexAttributeDescriptions),
      vsg::InputAssemblyState::create(),
      vsg::RasterizationState::create(),
      vsg::MultisampleState::create(),
      vsg::ColorBlendState::create(),
      vsg::DepthStencilState::create() };

    auto pipelineLayout = vsg::PipelineLayout::create(vsg::DescriptorSetLayouts{ descriptorSetLayout}, pushConstantRanges);
    auto graphicsPipeline = vsg::GraphicsPipeline::create(pipelineLayout, vsg::ShaderStages{ vertexShader, fragmentShader }, pipelineStates);
    auto bindGraphicsPipeline = vsg::BindGraphicsPipeline::create(graphicsPipeline);

    // create texture image and associated DescriptorSets and binding
    auto clampToEdge_sampler = vsg::Sampler::create();
    clampToEdge_sampler->addressModeU = VK_SAMPLER_ADDRESS_MODE_CLAMP_TO_EDGE;
    clampToEdge_sampler->addressModeV = VK_SAMPLER_ADDRESS_MODE_CLAMP_TO_EDGE;

    auto texture = vsg::DescriptorImage::create(clampToEdge_sampler, textureData, 0, 0, VK_DESCRIPTOR_TYPE_COMBINED_IMAGE_SAMPLER);

    auto descriptorSet = vsg::DescriptorSet::create(descriptorSetLayout, vsg::Descriptors{ texture });
    auto bindDescriptorSet = vsg::BindDescriptorSet::create(VK_PIPELINE_BIND_POINT_GRAPHICS, graphicsPipeline->layout, 0, descriptorSet);

    // create StateGroup as the root of the scene/command graph to hold the GraphicsProgram, and binding of Descriptors to decorate the whole graph
    auto scenegraph = vsg::StateGroup::create();
    scenegraph->add(bindGraphicsPipeline);
    scenegraph->add(bindDescriptorSet);

    auto geometry = get_volume_rendering();
    vsg::dvec3 position(0.0f, 0.0f, 0.0f);
    auto transform = vsg::MatrixTransform::create(vsg::translate(position));

    // add geometry
    transform->addChild(geometry);

    scenegraph->addChild(transform);

    // create the viewer and assign window(s) to it
    auto viewer = vsg::Viewer::create();

    auto window = vsg::Window::create(windowTraits);
    if (!window)
    {
        std::cout << "Could not create windows." << std::endl;
        return 1;
    }

    viewer->addWindow(window);

    // camera related details
    auto viewport = vsg::ViewportState::create(0, 0, window->extent2D().width, window->extent2D().height);
    auto perspective = vsg::Perspective::create(60.0, static_cast<double>(window->extent2D().width) / static_cast<double>(window->extent2D().height), 0.1, 10.0);
    auto lookAt = vsg::LookAt::create(vsg::dvec3(1.0, 1.0, 1.0), vsg::dvec3(0.0, 0.0, 0.0), vsg::dvec3(0.0, 0.0, 1.0));
    auto camera = vsg::Camera::create(perspective, lookAt, viewport);


    auto commandGraph = vsg::createCommandGraphForView(window, camera, scenegraph);
    viewer->assignRecordAndSubmitTaskAndPresentation({ commandGraph });

    // compile the Vulkan objects
    viewer->compile();

    // assign a CloseHandler to the Viewer to respond to pressing Escape or press the window close button
    viewer->addEventHandler(vsg::CloseHandler::create(viewer));
    viewer->addEventHandler(vsg::Trackball::create(camera));

    // main frame loop
    while (viewer->advanceToNextFrame())
    {
        // pass any events into EventHandlers assigned to the Viewer
        viewer->handleEvents();

        viewer->update();

        viewer->recordAndSubmit();

        viewer->present();
    }

    // clean up done automatically thanks to ref_ptr<>
    return 0;
}

[robertosfield]

This might a good example for vsgExamples, so please submit this for inclusion if you hare happy to do this. A shared example would help with debugging.

What results are you seeing right now?

One thing that might be an issue is the VSG by default uses reverse depth so the depth fragment go from 0 at the far plane to 1.0 at the near plane, this is opposite the OSG default.

[Joaopmoliveira]

This is the behavior that I am observing. It does seem that I am rendering only the surface of the volume itself. I already tried to recreate the accumulate functions in OSG and the behavior is still not what I expected.

I will try to read the thread you wrote in the google forum to understand the implications for the frag shader of reversing the depth.

[Joaopmoliveira]

I will try to put together an example in vsgExamples where the code can be refined. Thank you!

[robertosfield]

Thanks for the example. I have merged the PR as Joaopmoliveira-master branch:

/~https://github.com/vsg-dev/vsgExamples/tree/Joaopmoliveira-master

I have renamed the directory from osg_reproduction to vsgvolume with the hope that this can evolve into an official example.

I have also changed the code to by default check the VSG_FILE_PATH variable for volume.vert and volume.frag shaders and have placed these in the vsgExampels/data/shaders directory. This means we can edit the .vert and .frag shaders and just re-run the example.

I will now have a tinker to see what might be going on.

[robertosfield]

I have changed the face culling from the default back face culling to front face culling, this is one of the trick that the osgVolume uses to get the start point of the rate to march back towards the eye point and wasn't yet implemented in vsgvolume:

vsg-dev/vsgExamples@4f47b59

I now see the brick in the default camera position bit it doesn't yet do the ray tracing correctly. Looking at the vertex and fragments shaders I now don't think the reverse depth will be an issue and won't require any special handling as the ray is generated in objects space rather than in clip space.

[Joaopmoliveira]

I have two obvious problems with the shaders I implemented. The texgen homogeneous matrix is just the identity in my implementation (defined in the vertex shader) which is wrong. In OpenGL you used the gl_ObjectPlaneS ,T, ... to construct the matrix, but I could not figure out how to properly port this to VSG.

Another question I have is how to properly define the texture data, .i.e. should it be a float or a r32b32g32a32 array?

[Joaopmoliveira]

And the values

    float TransparencyValue = 0.4;
    float AlphaFuncValue = 0.7;
    float SampleDensityValue = 10;

should be uniform float value, defined and binded to the pipeline from the VSG code and not in the shaders obviously.

[robertosfield]

I have some half reasonable now checked in to the Joaopmoliveira-master branch:

I'm currently trying to fix an automated build issue, once that's done I'll merge these changes with vsgExamples master. From there there is other work that could be done.

I've added the ability to load a 3D texture image on the command line but it doesn't work yet because the vsgXchange DDS reader writer doesn't yet support 3D texture images, and we don't have a DIACOM loader like the OSG has. The DDS reader should be a first step to get working. I'm thinking a vsgVolume library that adds an DIACOM loader and other volume related features like osgVolume does would be an appropriate thing to do longer term.

[robertosfield]

Changes are now merged with vsgExamples master.

[Joaopmoliveira]

Good morning @robertosfield
I would like to continue with this line of work and implement a prototype of proper volume rendering, where currently, the example committed ignores:

1. z-depth of other objects rendered with VSG (depending on the ordering of the rendered objects, a plane with four corners cutting the volume in half completely cuts the region closest to the camera instead of being interleaved with the said plane)
2. transfer functions to inject color depending on the grayscale value

I have looked at the Openscene graph to implement the transfer function functionally, but for my life, I can't figure out how to deal with the z-depth problem. Solutions I have looked at:

Ordering could be enforced such that all other traditional geometries are first rendered, and only afterward are the volumes rendered. Do you think I should add a separate renderpass, similar to what you implemented with Im-GUI? Even if this extra render pass is implemented, the interleaving problem mentioned in 1) remains unsolved to the best of my knowledge.

I would like, if possible, to get your first impression on what I would have to do to have this feature, as well as how you would go about to solving it. From there I will try to get a working prototype with your insights so that we can add it to your examples, and hopefully to your amazing vsg library itself.

On another note, for the DICOM reader, would you happen to have a preferred library to do so? I currently use ITK for most of my medical applications, but I would understand if the size of the library makes it an undesirable feature to add to vsgXchange.

[robertosfield]

A long time back when I did volume rendering work with the OpenSceneGraph, to handle the rendering of 3d meshes in the same scene as a volume I have previously used a couple of approaches:

1. Render the meshes to a depth buffer, then render the 3d volume cubes front face and then read the depth buffer at that fragment position to get the depth to trace through the volume to.
2. Render the meshes with a volume rendering fragment shader that both lights the surface and does the ray trace of the volume, then render the back faces of the volume and ray trace to render the rest of the 3d volume outside of the fragments affected by the mesh.

The transfer function is simply a 1D texture on the GPU, to conveniently set that up I've previously using a std::map<double, vec4> to provide the input data and then interpolate the data to fill in the 1D texture.

For the OSG project I wrote a DICOM plugin using the DCMTK library to read the files. I haven't done any serious work on the volume rendering side for a decade so I don't know what would be the best solution now. The OSG's dicom plugin could serve as a reference for any new implementation. The basics of image conversion would be similar.

All this functionality could be provided as an vsgVolume library like has happened with vsgPoints, vsgImGui etc. However, it's not something I can personally embark on without funding for the effort put it.

[Joaopmoliveira]

The example would still be welcome in the project, correct? At least to add to vsgExamples for future inspiration. If not, I understand entirely.
Otherwise I can start the implementation, and any feedback would be helpful! I wrote a wrapper that does some of this functionality for me, so its not of the upmost importance to have this shipped with VSG itself.

[robertosfield]

Your are welcome to continue to extend the vsgvolume example, or even begin work on open source volume rendering library :-)