Further advances in nuklear port.

Rendering now looks similar to Qt version, needs a few tweaks:
* Proper multisampling
* Background

Missing features:
* Barcharts
* Interactivity (e.g. brushing/linking in all objects)
* History view of interactions

1 files changed, 429 insertions, 21 deletions

diff --git a/voronoisplat.cpp b/voronoisplat.cpp
index 3bbe7cf..1bcab69 100644
--- a/voronoisplat.cpp
+++ b/voronoisplat.cpp
@@ -4,12 +4,150 @@
 
 #include "colormap.h"
 #include "scatterplot.h"
+#include "shader.h"
 // #include "skelft.h"
 
 static const float DEFAULT_ALPHA = 5.0f;
 static const float DEFAULT_BETA = 20.0f;
 
-static const int SAMPLES = 128;
+static const char *programVoronoiVertexShader = R"EOF(#version 330
+// Single triangle strip quad generated entirely on the vertex shader.
+// Simply do glDrawArrays(GL_TRIANGLE_STRIP, 0, 4) and the shader
+// generates 4 points from gl_VertexID. No Vertex Attributes are
+// required.
+
+precision mediump float;
+
+uniform mat4 transform;
+
+layout(location = 0) in vec2 site_pos;
+// layout(location = 1) in vec4 site_color;
+
+out vec2 site;
+// out vec4 color;
+
+void main(void)
+{
+    vec2 uv;
+    uv.x = (gl_VertexID & 1);
+    uv.y = ((gl_VertexID >> 1) & 1);
+    gl_Position = vec4(uv * 2.0 - 1.0, 0.0, 1.0);
+    vec4 site4d = transform * vec4(site_pos.xy, 0.0, 1.0);
+    // 0.5 otherwise our maps are always scaled twice
+    site = 0.5 * site4d.xy;
+    // color = site_color;
+}
+)EOF";
+
+static const char *programVoronoiFragmentShader = R"EOF(#version 330
+
+precision mediump float;
+
+uniform vec2 resolution;
+uniform float rad_blur;
+uniform float rad_max;
+
+// in vec4 color;
+in vec2 site;
+out vec4 color;
+
+void main(void) {
+    float dt = length(gl_FragCoord.xy - site);
+    float maxlen = length(resolution);
+    float normalised_dt = dt / maxlen;
+    // float radius = rad_max + rad_blur;
+    // float normalised_radius = radius / maxlen;
+    gl_FragDepth = normalised_dt;
+    color = vec4(dt, 0.0, 0.0, 0.0);
+}
+)EOF";
+
+static const char *program1VertexShader = R"EOF(#version 330
+
+uniform float rad_blur;
+uniform float rad_max;
+uniform mat4 transform;
+
+in vec2 vert;
+in float scalar;
+
+out float value;
+
+void main() {
+    gl_PointSize = 2.0 * (rad_max + rad_blur);
+    gl_Position = transform * vec4(vert, 0.0, 1.0);
+    value = scalar;
+}
+)EOF";
+
+static const char *program1FragmentShader = R"EOF(#version 330
+
+uniform float rad_blur;
+uniform float rad_max;
+uniform sampler2D siteDT;
+
+in float value;
+
+layout (location = 0) out vec4 fragColor;
+
+void main() {
+    float dt = texelFetch(siteDT, ivec2(gl_FragCoord.xy), 0).r;
+    if (dt > rad_max) {
+        discard;
+    } else {
+        vec2 point = gl_PointCoord - vec2(0.5, 0.5);
+        float d2 = dot(point, point);
+        float radius = rad_max + rad_blur;
+        // float r = 2.0 * distance(gl_PointCoord, vec2(0.5, 0.5)) * radius;
+        // float r2 = r * r;
+        float r2 = 4.0 * d2 * radius * radius;
+        float dt_blur = dt + rad_blur;
+        float dt_blur2 = dt_blur * dt_blur;
+        if (r2 > dt_blur2) {
+            discard;
+        } else {
+            float w = exp(-5.0 * r2 / dt_blur2);
+            fragColor = vec4(w * value, w, 0.0, 0.0);
+        }
+    }
+    // fragColor = vec4(0.0, 0.0, 0.0, 0.0);
+}
+)EOF";
+
+static const char *program2VertexShader = R"EOF(#version 330
+
+in vec2 vert;
+
+void main() {
+    gl_Position = vec4(vert, 0.0, 1.0);
+}
+)EOF";
+
+static const char *program2FragmentShader = R"EOF(#version 330
+
+uniform sampler2D siteDT;
+uniform sampler2D accumTex;
+uniform sampler2D colormap;
+uniform float rad_max;
+
+layout (location = 0) out vec4 fragColor;
+
+vec3 getRGB(float value) {
+    return texture(colormap, vec2(mix(0.005, 0.995, value), 0)).rgb;
+}
+
+void main() {
+    float dt = texelFetch(siteDT, ivec2(gl_FragCoord.xy), 0).r;
+    if (dt > rad_max) {
+        discard;
+    } else {
+        vec4 accum = texelFetch(accumTex, ivec2(gl_FragCoord.xy), 0);
+        float value = accum.g > 0.0 ? accum.r / accum.g : 0.0;
+        // float value = (accum.g > 1.0) ? (accum.r - 1.0) / (accum.g - 1.0) : 0.0;
+        fragColor = vec4(getRGB(value), 1.0 - dt / rad_max);
+    }
+}
+)EOF";
 
 static int nextPow2(int n)
 {
@@ -29,11 +167,37 @@ VoronoiSplat::VoronoiSplat()
     , m_sitesChanged(false)
     , m_valuesChanged(false)
     , m_colorScaleChanged(false)
+    , m_redraw(false)
 {
+    std::fill(&m_transform[0][0], &m_transform[0][0] + 16, 0.0f);
+    m_transform[3][3] = 1.0f;
+
+    glGenFramebuffers(1, &m_voronoiFBO);
+    glGenFramebuffers(1, &m_preFBO);
+    glGenFramebuffers(1, &m_FBO);
+
+    setupShaders();
+    setupVAOs();
+    setupTextures();
+}
+
+VoronoiSplat::~VoronoiSplat()
+{
+    glDeleteFramebuffers(1, &m_FBO);
+    glDeleteFramebuffers(1, &m_preFBO);
+    glDeleteFramebuffers(1, &m_voronoiFBO);
 }
 
 void VoronoiSplat::update()
 {
+    m_redraw = true;
+
+    if (m_sitesChanged) {
+        updateSites();
+    }
+    if (m_valuesChanged) {
+        updateValues();
+    }
 }
 
 void VoronoiSplat::setSites(const arma::mat &points)
@@ -48,16 +212,11 @@ void VoronoiSplat::setSites(const arma::mat &points)
     }
 
     // Copy 'points' to internal data structure(s)
-    // Coords are packed into 'm_sites' as [ x1, y1, x2, y2, ... ]
-    m_sites.resize(2*points.n_rows);
-    const double *col = points.colptr(0);
-    for (unsigned i = 0; i < points.n_rows; i++) {
-        m_sites[2*i] = col[i];
-    }
-
-    col = points.colptr(1);
+    m_sites.resize(points.n_rows);
+    const double *col_x = points.colptr(0);
+    const double *col_y = points.colptr(1);
     for (unsigned i = 0; i < points.n_rows; i++) {
-        m_sites[2*i + 1] = col[i];
+        m_sites[i].set(col_x[i], col_y[i]);
     }
 
     setSitesChanged(true);
@@ -67,7 +226,7 @@ void VoronoiSplat::setSites(const arma::mat &points)
 void VoronoiSplat::setValues(const arma::vec &values)
 {
     if (values.n_elem == 0
-        || (m_sites.size() != 0 && values.n_elem != m_sites.size() / 2)) {
+        || (m_sites.size() != 0 && values.n_elem != m_sites.size())) {
         return;
     }
 
@@ -80,16 +239,6 @@ void VoronoiSplat::setValues(const arma::vec &values)
     update();
 }
 
-void VoronoiSplat::setColorScale(std::shared_ptr<const ColorScale> scale)
-{
-    m_cmap.resize(SAMPLES * 3);
-    scale->sample(SAMPLES, m_cmap.begin());
-    colorScaleChanged(scale);
-
-    setColorScaleChanged(true);
-    update();
-}
-
 void VoronoiSplat::setScale(const LinearScale<float> &sx,
                             const LinearScale<float> &sy)
 {
@@ -99,6 +248,13 @@ void VoronoiSplat::setScale(const LinearScale<float> &sx,
     update();
 }
 
+void VoronoiSplat::setColormap(GLuint texture)
+{
+    m_colormapTex = texture;
+    colormapChanged(texture);
+    update();
+}
+
 void VoronoiSplat::setAlpha(float alpha)
 {
     m_alpha = alpha;
@@ -113,6 +269,258 @@ void VoronoiSplat::setBeta(float beta)
     update();
 }
 
+void VoronoiSplat::setSize(size_t width, size_t height)
+{
+    m_width = width;
+    m_height = height;
+
+    resizeTextures();
+}
+
+void VoronoiSplat::setupShaders()
+{
+    m_program1 = std::make_unique<Shader>(
+        program1VertexShader,
+        program1FragmentShader);
+    m_program2 = std::make_unique<Shader>(
+        program2VertexShader,
+        program2FragmentShader);
+    m_programVoronoi = std::make_unique<Shader>(
+        programVoronoiVertexShader,
+        programVoronoiFragmentShader);
+}
+
+void VoronoiSplat::setupVAOs()
+{
+    // sitesVAO: VBOs 0 & 1 are for sites & their values (init'd later)
+    glGenBuffers(3, m_VBOs);
+    glGenVertexArrays(1, &m_sitesVAO);
+    glBindVertexArray(m_sitesVAO);
+    glBindBuffer(GL_ARRAY_BUFFER, m_VBOs[0]);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, nullptr);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_VBOs[1]);
+    glEnableVertexAttribArray(1);
+    glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 0, nullptr);
+    glBindVertexArray(0);
+
+    // 2ndPassVAO: VBO 2 is a quad mapping the final texture to the framebuffer
+    glGenVertexArrays(1, &m_2ndPassVAO);
+    glBindVertexArray(m_2ndPassVAO);
+    GLfloat verts[] = { -1.0f, -1.0f, -1.0f,  1.0f,
+                         1.0f, -1.0f,  1.0f,  1.0f };
+    glBindBuffer(GL_ARRAY_BUFFER, m_VBOs[2]);
+    glBufferData(GL_ARRAY_BUFFER, sizeof(verts), verts, GL_STATIC_DRAW);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, nullptr);
+    glBindVertexArray(0);
+
+    glGenVertexArrays(1, &m_voronoiVAO);
+    glBindVertexArray(m_voronoiVAO);
+    glGenBuffers(1, &m_dtVBO);
+    glBindBuffer(GL_ARRAY_BUFFER, m_dtVBO);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, nullptr);
+    glVertexAttribDivisor(0, 1);
+    glBindVertexArray(0);
+}
+
+void VoronoiSplat::setupTextures()
+{
+    glGenTextures(2, m_textures);
+
+    // Texture where output is drawn to
+    glGenTextures(1, &m_outTex);
+
+    // Voronoi diagram is built by relying on depth tests on GPU
+    glGenTextures(1, &m_voronoiDepthTex);
+}
+
+void VoronoiSplat::resizeTextures()
+{
+    // textures[0] stores the DT values for each pixel
+    glBindTexture(GL_TEXTURE_2D, m_textures[0]);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, m_width, m_height, 0, GL_RED,
+                 GL_FLOAT, 0);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+
+    // textures[1] is the result of the first pass
+    glBindTexture(GL_TEXTURE_2D, m_textures[1]);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, m_width, m_height, 0, GL_RGBA,
+                 GL_FLOAT, 0);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glBindTexture(GL_TEXTURE_2D, 0);
+
+    glBindTexture(GL_TEXTURE_2D, m_voronoiDepthTex);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT, m_width, m_height, 0, GL_DEPTH_COMPONENT, GL_UNSIGNED_BYTE, NULL);
+    glBindTexture(GL_TEXTURE_2D, 0);
+
+    glBindTexture(GL_TEXTURE_2D, m_outTex);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, m_width, m_height, 0, GL_RGBA,
+                 GL_FLOAT, 0);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
+    glBindTexture(GL_TEXTURE_2D, 0);
+}
+
+void VoronoiSplat::updateSites()
+{
+    glBindBuffer(GL_ARRAY_BUFFER, m_VBOs[0]);
+    glBufferData(GL_ARRAY_BUFFER, m_sites.size() * sizeof(vec2),
+                 m_sites.data(), GL_DYNAMIC_DRAW);
+
+    // Compute DT values for the new positions
+    // computeDT();
+    float padding = Scatterplot::PADDING;
+    m_sx.setRange(static_cast<float>(padding),
+                  static_cast<float>(m_width - padding));
+    m_sy.setRange(static_cast<float>(m_height - padding),
+                  static_cast<float>(padding));
+    updateTransform4x4(m_sx, m_sy, m_transform);
+
+    GLint originalFBO;
+    glGetIntegerv(GL_FRAMEBUFFER_BINDING, &originalFBO);
+    glBindFramebuffer(GL_FRAMEBUFFER, m_voronoiFBO);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D,
+                           m_textures[0], 0);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_TEXTURE_2D,
+                           m_voronoiDepthTex, 0);
+
+    glViewport(0, 0, m_width, m_height);
+    m_programVoronoi->use();
+    m_programVoronoi->setUniform("transform", m_transform);
+    m_programVoronoi->setUniform("rad_max", m_beta);
+    m_programVoronoi->setUniform("rad_blur", m_alpha);
+    GLfloat resolution[] = {
+        static_cast<GLfloat>(m_width),
+        static_cast<GLfloat>(m_height),
+    };
+    m_programVoronoi->setUniform2dArray("resolution", resolution, 1);
+    glEnable(GL_DEPTH_TEST);
+    glEnable(GL_BLEND);
+    // glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ZERO);
+    glBlendFunc(GL_ONE, GL_ZERO);
+
+    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+    glBindVertexArray(m_voronoiVAO);
+    glBindBuffer(GL_ARRAY_BUFFER, m_dtVBO);
+    glBufferData(GL_ARRAY_BUFFER, m_sites.size() * sizeof(vec2),
+                 m_sites.data(), GL_DYNAMIC_DRAW);
+    glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, m_sites.size());
+    glBindVertexArray(0);
+    m_programVoronoi->release();
+    glDisable(GL_DEPTH_TEST);
+    glBindFramebuffer(GL_FRAMEBUFFER, originalFBO);
+
+    // Update transform used when drawing sites
+    updateTransform();
+
+    m_sitesChanged = false;
+}
+
+void VoronoiSplat::updateValues()
+{
+    glBindBuffer(GL_ARRAY_BUFFER, m_VBOs[1]);
+    glBufferData(GL_ARRAY_BUFFER, m_values.size() * sizeof(GLfloat),
+                 m_values.data(), GL_DYNAMIC_DRAW);
+
+    m_valuesChanged = false;
+}
+
+void VoronoiSplat::updateTransform()
+{
+    float padding = Scatterplot::PADDING;
+    float offsetX = padding / m_width, offsetY = padding / m_height;
+    updateTransform4x4(m_sx, m_sy, offsetX, offsetY, m_transform);
+}
+
+void VoronoiSplat::draw()
+{
+    if (!m_redraw) {
+        return;
+    }
+
+    int originalFBO;
+    glGetIntegerv(GL_FRAMEBUFFER_BINDING, &originalFBO);
+
+    // First, we draw to an intermediate texture, which is used as input for the
+    // second pass
+    glBindFramebuffer(GL_FRAMEBUFFER, m_preFBO);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
+                           GL_TEXTURE_2D, m_textures[1], 0);
+
+    glViewport(0, 0, m_width, m_height);
+    m_program1->use();
+    m_program1->setUniform("rad_max", m_beta);
+    m_program1->setUniform("rad_blur", m_alpha);
+    m_program1->setUniform("transform", m_transform);
+
+    glActiveTexture(GL_TEXTURE0);
+    glBindTexture(GL_TEXTURE_2D, m_textures[0]);
+    m_program1->setUniform("siteDT", 0);
+
+    // glEnable(GL_POINT_SPRITE);
+    glEnable(GL_PROGRAM_POINT_SIZE);
+    glEnable(GL_BLEND);
+    glBlendFunc(GL_ONE, GL_ZERO);
+
+    glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+
+    glBindVertexArray(m_sitesVAO);
+    glDrawArrays(GL_POINTS, 0, m_sites.size());
+    glBindVertexArray(0);
+    m_program1->release();
+    glDisable(GL_PROGRAM_POINT_SIZE);
+
+    // Second pass
+    m_program2->use();
+    m_program2->setUniform("rad_max", m_beta);
+
+    glActiveTexture(GL_TEXTURE0);
+    glBindTexture(GL_TEXTURE_2D, m_textures[0]);
+    m_program2->setUniform("siteDT", 0);
+    glActiveTexture(GL_TEXTURE1);
+    glBindTexture(GL_TEXTURE_2D, m_textures[1]);
+    m_program2->setUniform("accumTex", 1);
+    glActiveTexture(GL_TEXTURE2);
+    glBindTexture(GL_TEXTURE_2D, m_colormapTex);
+    m_program2->setUniform("colormap", 2);
+
+    glBindFramebuffer(GL_FRAMEBUFFER, m_FBO);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
+                           GL_TEXTURE_2D, m_outTex, 0);
+    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+    // TODO: know beforehand which color to be used as transparent for
+    // blending. We currently assume we always plot to a white
+    // background.
+    glClearColor(1.0f, 1.0f, 1.0f, 0.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+
+    glBindVertexArray(m_2ndPassVAO);
+    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
+    glBindVertexArray(0);
+    m_program2->release();
+
+    /*
+    GLenum status = glCheckFramebufferStatus(GL_FRAMEBUFFER);
+    if (status != GL_FRAMEBUFFER_COMPLETE) {
+        checkGLError("draw():");
+    }
+    */
+
+    glBindFramebuffer(GL_FRAMEBUFFER, originalFBO);
+    m_redraw = false;
+}
+
 // ----------------------------------------------------------------------------
 /*
 class VoronoiSplatRenderer