volume_canvas.cpp

Go to the documentation of this file.

#include "volume_canvas.h"

#include <wx/dcclient.h>
#include <wx/image.h>

#include <fstream>
#include <vector>
#include <algorithm>
#include <stdexcept>
#include <utility>

BEGIN_EVENT_TABLE(VolumeCanvas, wxGLCanvas)
    EVT_PAINT(VolumeCanvas::OnPaint)
    EVT_SIZE(VolumeCanvas::OnSize)
    EVT_ERASE_BACKGROUND(VolumeCanvas::OnEraseBackground)
END_EVENT_TABLE()

namespace
{
    std::pair<float, float> ComputeScaleFactors(float a1, float a2)
    {
        const float x = a2 / a1;
        if (x > 1.0f) return std::make_pair(1.0f, 1.0f / x);
        else return std::make_pair(x, 1.0f);
    }

    int attrib_list[] = { WX_GL_RGBA, WX_GL_DOUBLEBUFFER, 0 };
}

VolumeCanvas::VolumeCanvas(wxWindow * parent, wxWindowID id)
    : wxGLCanvas(parent, id, wxDefaultPosition, wxDefaultSize, wxFULL_REPAINT_ON_RESIZE, wxGLCanvasName,attrib_list)
    , m_volume_texture(0)
    , m_transfer_texture(0)
    , m_glstatus(GLS_INIT_PENDING)
    , m_volprog(0 )
    , m_viewmode(VM_SLICE_X)
    , m_sizex(0)
    , m_sizey(0)
    , m_sizez(0)
    , m_slice(0)
    , m_rotx(0)
    , m_roty(0)
    , m_rotz(0)
{
}

void VolumeCanvas::UpdateTransferFunction(TransferFunctionSampler *sampler)
{
    std::vector<float> data;

    // make sure density 0 is mapped to a completely transparent color
    data.push_back(0.0f);
    data.push_back(0.0f);
    data.push_back(0.0f);
    data.push_back(0.0f);

    for (int i = 1; i < 4096; ++i) {
        const TransferFunctionSampler::Color c = sampler->MapDensity(i / 4095.0f);
        data.push_back(c.r);
        data.push_back(c.g);
        data.push_back(c.b);
        data.push_back(c.a);
    }

    glBindTexture(GL_TEXTURE_1D, m_transfer_texture);
    glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA, 4096, 0, GL_RGBA, GL_FLOAT, &data[0]);

    Refresh();
}

void VolumeCanvas::Render3dView() const
{
    // cube vertex data
    static const GLfloat vdata[] =
    {
         1.0f,  1.0f, -1.0f,
         1.0f, -1.0f, -1.0f,
        -1.0f, -1.0f, -1.0f,
        -1.0f,  1.0f, -1.0f,
         1.0f,  1.0f,  1.0f,
         1.0f, -1.0f,  1.0f,
        -1.0f, -1.0f,  1.0f,
        -1.0f,  1.0f,  1.0f
    };

    // cube index data
    static const GLushort idata[] =
    {
        0, 1, 2, 3,
        4, 7, 6, 5,
        0, 4, 5, 1,
        1, 5, 6, 2,
        2, 6, 7, 3,
        4, 0, 3, 7
    };

    const double aspect = static_cast<double>(m_viewport_width) / m_viewport_height;

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    // TODO: let the user specify fovy and possibly other parameters?
    gluPerspective(54.0, aspect, 0.1, 100.0);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    // move the camera a bit away
    glTranslatef(0.0f, 0.0f, -3.0f);

    // orient the volume according to the settings
    glRotatef(m_rotx, 1.0f, 0.0f, 0.0f);
    glRotatef(m_roty, 0.0f, 1.0f, 0.0f);
    glRotatef(m_rotz, 0.0f, 0.0f, 1.0f);

    // compute scaling factor for the box so that the aspect ratio is correct
    const float ms = std::max(std::max(GetSizeX(), GetSizeY()), GetSizeZ());
    const float sx = GetSizeX() / ms;
    const float sy = GetSizeY() / ms;
    const float sz = GetSizeZ() / ms;
    glScalef(sx, sy, sz);

    m_volprog->SetMode(VolProgram::M_3D_COLOR_CUBE);
    m_volprog->Use();

    GLuint tex[2];
    glGenTextures(2, tex);

    // draw the cube

    glEnable(GL_CULL_FACE);
    glEnableClientState(GL_VERTEX_ARRAY);
    glVertexPointer(3, GL_FLOAT, 0, vdata);

    // draw the back faces
    glCullFace(GL_FRONT);
    glDrawElements(GL_QUADS, 24, GL_UNSIGNED_SHORT, idata);
    glBindTexture(GL_TEXTURE_2D, tex[0]);
    glCopyTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 0, 0, m_viewport_width, m_viewport_height, 0);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

    // draw the front faces
    glCullFace(GL_BACK);
    glDrawElements(GL_QUADS, 24, GL_UNSIGNED_SHORT, idata);
    glBindTexture(GL_TEXTURE_2D, tex[1]);
    glCopyTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, 0, 0, m_viewport_width, m_viewport_height, 0);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);

    glDisableClientState(GL_VERTEX_ARRAY);
    glDisable(GL_CULL_FACE);

    m_volprog->SetMode(VolProgram::M_3D_RAYCAST);
    m_volprog->SetBackAndFrontTextures(tex[0], tex[1]);
    m_volprog->SetVolumeTexture(m_volume_texture);
    m_volprog->SetTransferFunctionTexture(m_transfer_texture);
    m_volprog->Use();

    // draw a fullscreen quad

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glBegin(GL_QUADS);
    glVertex2f(-1.0f, -1.0f);
    glVertex2f( 1.0f, -1.0f);
    glVertex2f( 1.0f,  1.0f);
    glVertex2f(-1.0f,  1.0f);
    glEnd();

    m_volprog->SetBackAndFrontTextures(0, 0);

    glDeleteTextures(2, tex);
}

void VolumeCanvas::RenderSliceView() const
{
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    const float aspect = static_cast<float> (m_viewport_width) / m_viewport_height;
    float slice;
    std::pair<float, float> scale;

    m_volprog->SetMode(VolProgram::M_SLICE);
    m_volprog->SetVolumeTexture(m_volume_texture);
    m_volprog->SetTransferFunctionTexture(m_transfer_texture);
    m_volprog->Use();

    switch (GetViewmode()) {
    case VM_SLICE_X:
        slice = static_cast<float> (GetSlice()) / GetSizeX();
        scale = ComputeScaleFactors(aspect, static_cast<float> (GetSizeY()) / GetSizeZ());
        glScalef(scale.first, scale.second, 1.0f);
        glBegin(GL_QUADS);
        glTexCoord3f(slice, 0.0f, 0.0f);
        glVertex2f(-1.0f, -1.0f);
        glTexCoord3f(slice, 1.0f, 0.0f);
        glVertex2f(1.0f, -1.0f);
        glTexCoord3f(slice, 1.0f, 1.0f);
        glVertex2f(1.0f, 1.0f);
        glTexCoord3f(slice, 0.0f, 1.0f);
        glVertex2f(-1.0f, 1.0f);
        glEnd();
        break;
    case VM_SLICE_Y:
        slice = static_cast<float> (GetSlice()) / GetSizeY();
        scale = ComputeScaleFactors(aspect, static_cast<float> (GetSizeX()) / GetSizeZ());
        glScalef(scale.first, scale.second, 1.0f);
        glBegin(GL_QUADS);
        glTexCoord3f(0.0f, slice, 0.0f);
        glVertex2f(-1.0f, -1.0f);
        glTexCoord3f(1.0f, slice, 0.0f);
        glVertex2f(1.0f, -1.0f);
        glTexCoord3f(1.0f, slice, 1.0f);
        glVertex2f(1.0f, 1.0f);
        glTexCoord3f(0.0f, slice, 1.0f);
        glVertex2f(-1.0f, 1.0f);
        glEnd();
        break;
    case VM_SLICE_Z:
        slice = static_cast<float> (GetSlice()) / GetSizeZ();
        scale = ComputeScaleFactors(aspect, static_cast<float> (GetSizeX()) / GetSizeY());
        glScalef(scale.first, scale.second, 1.0f);
        glBegin(GL_QUADS);
        glTexCoord3f(0.0f, 0.0f, slice);
        glVertex2f(-1.0f, -1.0f);
        glTexCoord3f(1.0f, 0.0f, slice);
        glVertex2f(1.0f, -1.0f);
        glTexCoord3f(1.0f, 1.0f, slice);
        glVertex2f(1.0f, 1.0f);
        glTexCoord3f(0.0f, 1.0f, slice);
        glVertex2f(-1.0f, 1.0f);
        glEnd();
        break;
    }
}

void VolumeCanvas::RenderToFile(const wxString& file_name, int width, int height)
{
    //TODO: Clean up RenderToFile ... it's ugly (though it works ... hopefully)

    if (!CheckGL()) return;

    GLuint id_rbf;
    GLuint id_fbo;

    //Create and bind Frame Buffer Object
    glGenFramebuffersEXT(1, &id_fbo);
    glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, id_fbo);

    //Create and bind Render Buffer
    glGenRenderbuffersEXT(1, &id_rbf);
    glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT, GL_RGBA, width, height);
    glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, id_rbf);

    //Attach renderbuffer to framebuffer
    glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_RENDERBUFFER_EXT, id_rbf);

    //If framebuffer status ok -> render offscreen image and save to file
    if (GL_FRAMEBUFFER_COMPLETE_EXT == glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT))
    {
        m_viewport_width = width;
        m_viewport_height = height;

        Render();

        unsigned char *rgb = (unsigned char*)malloc(width * height * 3);
        unsigned char *alpha = (unsigned char*)malloc(width * height);

        glPixelStorei(GL_PACK_ALIGNMENT, 1);
        glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, rgb);
        glReadPixels(0, 0, width, height, GL_ALPHA, GL_UNSIGNED_BYTE, alpha);
        glPixelStorei(GL_PACK_ALIGNMENT, 4);

        wxImage img(width, height);
        img.SetData(rgb);
        img.SetAlpha(alpha);
        img = img.Mirror(false);

        if (!img.SaveFile(file_name)) {
            // TODO: react on failure
        }
    } else {
        //TODO: react on frame/renderbuffer error
    }

    //Unbind frame buffer
    glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);

    //Delete renderbuffer
    glDeleteRenderbuffersEXT(1, &id_rbf);

    //Delete framebuffer
    glDeleteFramebuffersEXT(1, &id_fbo);
}

void VolumeCanvas::Render() const
{
    glViewport(0, 0, m_viewport_width, m_viewport_height);
    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    if (!m_volume_texture) return;

    if (GetViewmode() == VM_3D) {
        Render3dView();
    } else {
        RenderSliceView();
    }
}

bool VolumeCanvas::InitGL()
{
    if (glewInit() != GLEW_NO_ERROR) return false;
    if (!GLEW_VERSION_2_0) return false;

    try {
        m_volprog = new VolProgram("data/volprog");
    } catch (...) {
        return false;
    }

    glGenTextures(1, &m_transfer_texture);
    glBindTexture(GL_TEXTURE_1D, m_transfer_texture);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);

    return true;
}

bool VolumeCanvas::CheckGL()
{
#ifndef __WXMOTIF__
    if (!GetContext()) return false;
#endif

    SetCurrent();
    if (m_glstatus == GLS_INIT_PENDING) m_glstatus = InitGL() ? GLS_OK : GLS_INIT_FAILED;
    if (m_glstatus != GLS_OK) return false;
    return true;
}

void VolumeCanvas::OnPaint(wxPaintEvent & event)
{
    wxPaintDC dc(this);
    if (!CheckGL()) return;

    // client size is used in the Render() method
    GetClientSize(&m_viewport_width, &m_viewport_height);
    Render();
    SwapBuffers();
}

void VolumeCanvas::OnSize(wxSizeEvent &event)
{
    wxGLCanvas::OnSize(event);
    Update();
}

void VolumeCanvas::OnEraseBackground(wxEraseEvent& event)
{
    // Do nothing, to avoid flashing.
}

namespace
{
    struct VolumeData {
        unsigned short grad_x;
        unsigned short grad_y;
        unsigned short grad_z;
        unsigned short value;
    };
}

void VolumeCanvas::OpenFile(const char *filename)
{
    if (!CheckGL()) throw std::runtime_error("failed to initialise OpenGL");

    std::ifstream in(filename, std::ios::in | std::ios::binary);
    in.exceptions(std::ios::failbit | std::ios::badbit);

    unsigned short size_x, size_y, size_z;

    // load the header which contains size data
    in.read(reinterpret_cast<char*> (&size_x), sizeof(size_x));
    in.read(reinterpret_cast<char*> (&size_y), sizeof(size_y));
    in.read(reinterpret_cast<char*> (&size_z), sizeof(size_z));

    // correct the read data in case we are running on a big endian architecture
    wxUINT16_SWAP_ON_BE(size_x);
    wxUINT16_SWAP_ON_BE(size_y);
    wxUINT16_SWAP_ON_BE(size_z);

    GLint max_size;
    glGetIntegerv(GL_MAX_3D_TEXTURE_SIZE, &max_size);

    if (size_x < 16 || size_y < 16 || size_z < 16 || size_x > max_size || size_y > max_size || size_z > max_size) {
        throw std::runtime_error("dataset size exceeds hardware limits");
    }

    const unsigned count = size_x * size_y * size_z;
    std::vector<VolumeData> volume_data(count);

    std::auto_ptr<HistogramData> histogram(new HistogramData());

    // load the data
    for (unsigned int z = 0; z < size_z; ++z) {
        for (unsigned int y = 0; y < size_y; ++y) {
            for (unsigned int x = 0; x < size_x; ++x) {
                unsigned short value;
                in.read(reinterpret_cast<char*> (&value), sizeof(value));
                wxUINT16_SWAP_ON_BE(value);

                // the data in the data files only uses 12 bits.
                // extend it to use the full 16 bit range
                value <<= 4;

                volume_data[x + y * size_x + z * size_x * size_y].value = value;

                histogram->AddData(value);
            }
        }
    }

    // Helper class for easy 3D array access into a simple 1D vector
    struct ArrayAccess3d {
        std::vector<VolumeData> &data;
        unsigned short sx, sy, sz;

        ArrayAccess3d(std::vector<VolumeData> &data, unsigned short sx, unsigned short sy, unsigned short sz) :
            data(data), sx(sx), sy(sy), sz(sz)
        {
        }

        VolumeData& operator ()(int x, int y, int z)
        {
            // clamp the input coordinates to the border
            x = std::min(std::max(x, 0), sx - 1);
            y = std::min(std::max(y, 0), sy - 1);
            z = std::min(std::max(z, 0), sz - 1);

            // access the value in the 3D array
            return data[x + y * sx + z * sx * sy];
        }
    } volume(volume_data, size_x, size_y, size_z);

    // compute the gradient of the data
    for (unsigned short z = 0; z < size_z; ++z) {
        for (unsigned short y = 0; y < size_y; ++y) {
            for (unsigned short x = 0; x < size_x; ++x) {
                volume(x, y, z).grad_x = 0.5 * (volume(x + 1, y, z).value - volume(x - 1, y, z).value);
                volume(x, y, z).grad_y = 0.5 * (volume(x, y + 1, z).value - volume(x, y - 1, z).value);
                volume(x, y, z).grad_z = 0.5 * (volume(x, y, z + 1).value - volume(x, y, z - 1).value);
            }
        }
    }

    // clear any OpenGL errors
    while (glGetError() != GL_NO_ERROR)
        ;

    // TODO: handle the case where the 3D hardware does not natively support NPOT textures.
    // TODO: handle case where GL_RGBA16 does not actually have 16 bits per channel on current hardware.

    GLuint id = 0;
    glGenTextures(1, &id);
    glBindTexture(GL_TEXTURE_3D, id);
    glTexImage3D(GL_TEXTURE_3D, 0, GL_RGBA16, size_x, size_y, size_z, 0, GL_RGBA, GL_UNSIGNED_SHORT,
        &volume_data[0]);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
    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);

    if (glGetError() != GL_NO_ERROR) {
        glDeleteTextures(1, &id);
        throw std::runtime_error("OpenGL error");
    }

    // replace old histogram data
    m_histogramdata = histogram;

    // replace old texture by the new one
    glDeleteTextures(1, &m_volume_texture);
    m_volume_texture = id;

    m_sizex = size_x;
    m_sizey = size_y;
    m_sizez = size_z;

    // make sure the newly loaded dataset will be displayed
    Refresh();
}

---