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

19 files changed, 1659 insertions, 224 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..1bb097f
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*~
+*.o
\ No newline at end of file
diff --git a/Makefile b/Makefile
index b409bfa..058675a 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
-CFLAGS=-O2 -Iinclude `pkg-config --cflags armadillo glfw3`
+CFLAGS=-O0 -g -Iinclude `pkg-config --cflags armadillo glfw3`
 LIBS=-lGL `pkg-config --libs armadillo glfw3`
-OBJS=barchart.o brushinghandler.o colormap.o colorscale.o continuouscolorscale.o divergentcolorscale.o dist.o forcescheme.o geometry.o lamp.o manipulationhandler.o mapscalehandler.o measures.o projectionhistory.o quadtree.o scatterplot.o voronoisplat.o
+OBJS=barchart.o brushinghandler.o colormap.o colorscale.o continuouscolorscale.o divergentcolorscale.o dist.o forcescheme.o geometry.o lamp.o manipulationhandler.o mapscalehandler.o measures.o projectionhistory.o quadtree.o scale.o scatterplot.o shader.o voronoisplat.o
 
 all: pm
 
@@ -68,11 +68,19 @@ quadtree.o: quadtree.cpp quadtree.h geometry.h
 	@echo CC $@
 	@g++ $(CFLAGS) $(LIBS) -c quadtree.cpp
 
-scatterplot.o: scatterplot.cpp scatterplot.h colorscale.h scale.h continuouscolorscale.h geometry.h
+scale.o: scale.cpp scale.h
+	@echo CC $@
+	@g++ $(CFLAGS) $(LIBS) -c scale.cpp
+
+scatterplot.o: scatterplot.cpp scatterplot.h colorscale.h continuouscolorscale.h scale.h geometry.h
 	@echo CC $@
 	@g++ $(CFLAGS) $(LIBS) -c scatterplot.cpp
 
-voronoisplat.o: voronoisplat.cpp voronoisplat.h colormap.h scatterplot.h
+shader.o: shader.cpp shader.h
+	@echo CC $@
+	@g++ $(CFLAGS) $(LIBS) -c shader.cpp
+
+voronoisplat.o: voronoisplat.cpp voronoisplat.h colormap.h scatterplot.h shader.h
 	@echo CC $@
 	@g++ $(CFLAGS) $(LIBS) -c voronoisplat.cpp
 
diff --git a/colormap.cpp b/colormap.cpp
index 08a3d61..0a9b4cc 100644
--- a/colormap.cpp
+++ b/colormap.cpp
@@ -3,69 +3,38 @@
 #include <algorithm>
 #include <glad/gl.h>
 
-class ColormapTexture
-{
-public:
-    ColormapTexture(const std::vector<float> &cmap,
-                    Colormap::Orientation orientation = Colormap::Horizontal);
-    ~ColormapTexture();
-
-    bool hasAlphaChannel() const { return false; }
-    bool hasMipmaps() const { return false; }
-    void bind();
-    bool updateTexture();
-
-    int textureId() const { return m_texture; }
-    size_t width() const { return m_width; }
-    size_t height() const { return m_height; }
-
-    void setOrientation(Colormap::Orientation orientation);
-
-private:
-    Colormap::Orientation m_orientation;
-    // QSize m_size;
-    size_t m_width, m_height;
-    GLuint m_texture;
-    const std::vector<float> &m_cmap;
-};
-
-ColormapTexture::ColormapTexture(const std::vector<float> &cmap,
-                                 Colormap::Orientation orientation)
-    : m_cmap(cmap)
-    , m_width(0)
+Colormap::Colormap()
+    : m_width(0)
     , m_height(0)
+    , m_shouldUpdateTexture(false)
+    , m_orientation(Colormap::Horizontal)
 {
-    // Setup OpenGL texture
     glGenTextures(1, &m_texture);
     glBindTexture(GL_TEXTURE_2D, m_texture);
     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
     glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
-
-    setOrientation(orientation);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
 }
 
-ColormapTexture::~ColormapTexture()
+Colormap::~Colormap()
 {
     glDeleteTextures(1, &m_texture);
 }
 
-void ColormapTexture::bind()
+void Colormap::update()
 {
-    glBindTexture(GL_TEXTURE_2D, m_texture);
+    if (m_shouldUpdateTexture) {
+        updateTexture();
+    }
 }
 
-void ColormapTexture::setOrientation(Colormap::Orientation orientation)
+void Colormap::updateTexture()
 {
-    if (m_orientation == orientation) {
+    if (m_cmap.empty()) {
         return;
     }
 
-    m_orientation = orientation;
-    updateTexture();
-}
-
-bool ColormapTexture::updateTexture()
-{
     switch (m_orientation) {
     case Colormap::Horizontal:
         m_width = m_cmap.size() / 3;
@@ -78,23 +47,8 @@ bool ColormapTexture::updateTexture()
     }
 
     glBindTexture(GL_TEXTURE_2D, m_texture);
-    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, m_width, m_height,
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F, m_width, m_height,
             0, GL_RGB, GL_FLOAT, m_cmap.data());
-    return true;
-}
-
-Colormap::Colormap()
-    : m_shouldUpdateTexture(false)
-    , m_orientation(Colormap::Horizontal)
-{
-}
-
-Colormap::~Colormap()
-{
-}
-
-void Colormap::update()
-{
 }
 
 static void reverseCMap(std::vector<float> &cmap)
@@ -110,9 +64,9 @@ static void reverseCMap(std::vector<float> &cmap)
     }
 }
 
-void Colormap::setOrientation(Colormap::Orientation orientation)
+void Colormap::setOrientation(Colormap::Orientation _orientation)
 {
-    if (m_orientation == orientation) {
+    if (m_orientation == _orientation) {
         return;
     }
 
@@ -120,8 +74,8 @@ void Colormap::setOrientation(Colormap::Orientation orientation)
         reverseCMap(m_cmap);
     }
 
-    m_orientation = orientation;
-    m_shouldUpdateOrientation = true;
+    m_orientation = _orientation;
+    m_shouldUpdateTexture = true;
     update();
 }
 
@@ -138,38 +92,3 @@ void Colormap::setColorScale(std::shared_ptr<const ColorScale> scale)
     m_shouldUpdateTexture = true;
     update();
 }
-
-/*
-QSGNode *Colormap::newSceneGraph()
-{
-    QSGSimpleTextureNode *node = new QSGSimpleTextureNode;
-    m_texture = new ColormapTexture(m_cmap);
-
-    node->setTexture(m_texture);
-    node->setOwnsTexture(false);
-    node->setSourceRect(0, 0, m_texture->width(), m_texture->height());
-
-    return node;
-}
-
-QSGNode *Colormap::updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *)
-{
-    QSGNode *root = oldNode ? oldNode : newSceneGraph();
-
-    QSGSimpleTextureNode *node = static_cast<QSGSimpleTextureNode *>(root);
-    node->setRect(x(), y(), width(), height());
-
-    ColormapTexture *texture = static_cast<ColormapTexture *>(m_texture);
-    if (m_shouldUpdateOrientation) {
-        texture->setOrientation(m_orientation);
-        m_shouldUpdateOrientation = false;
-    }
-
-    if (m_shouldUpdateTexture) {
-        texture->updateTexture();
-        m_shouldUpdateTexture = false;
-    }
-
-    return root;
-}
-*/
diff --git a/colormap.h b/colormap.h
index 76577f5..7d26a86 100644
--- a/colormap.h
+++ b/colormap.h
@@ -5,6 +5,7 @@
 #include <vector>
 
 #include <armadillo>
+#include <glad/gl.h>
 #include <nod.hpp>
 
 #include "colorscale.h"
@@ -22,6 +23,10 @@ public:
     Colormap();
     ~Colormap();
 
+    GLuint texture() const { return m_texture; }
+    size_t width() const { return m_width; }
+    size_t height() const { return m_height; }
+
     void update();
     void setOrientation(Orientation orientation);
     Orientation orientation() const { return m_orientation; }
@@ -31,13 +36,12 @@ public:
 
     void setColorScale(std::shared_ptr<const ColorScale> scale);
 
-    // QSGNode *updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *);
-
 private:
-    // QSGNode *newSceneGraph();
+    void updateTexture();
 
-    // QSGDynamicTexture *m_texture;
-    bool m_shouldUpdateTexture, m_shouldUpdateOrientation;
+    size_t m_width, m_height;
+    GLuint m_texture;
+    bool m_shouldUpdateTexture;
     Orientation m_orientation;
     std::vector<float> m_cmap;
 };
diff --git a/colorscale.cpp b/colorscale.cpp
index 32d0f25..ec52d2d 100644
--- a/colorscale.cpp
+++ b/colorscale.cpp
@@ -37,15 +37,15 @@ Color::Color(float r, float g, float b, float a)
 
 void Color::getRgbF(float *r, float *g, float *b) const
 {
-    *r = round(static_cast<float>(this->r) / 255.0f);
-    *g = round(static_cast<float>(this->g) / 255.0f);
-    *b = round(static_cast<float>(this->b) / 255.0f);
+    *r = static_cast<float>(this->r) / 255.0f;
+    *g = static_cast<float>(this->g) / 255.0f;
+    *b = static_cast<float>(this->b) / 255.0f;
 }
 
 void Color::getRgbF(float *r, float *g, float *b, float *a) const
 {
     getRgbF(r, g, b);
-    *a = round(static_cast<float>(this->a) / 255.0f);
+    *a = static_cast<float>(this->a) / 255.0f;
 }
 
 void Color::setRgb(int r, int g, int b)
diff --git a/colorscale.h b/colorscale.h
index c50deb0..e2854bd 100644
--- a/colorscale.h
+++ b/colorscale.h
@@ -44,7 +44,7 @@ public:
     int numColors() const { return m_colors.size(); }
 
     template<typename OutputIterator>
-    void sample(int samples, OutputIterator it) const;
+    void sample(std::size_t samples, OutputIterator it) const;
 
     static Color lerp(const Color &c1, const Color &c2, float _t);
 
@@ -54,7 +54,7 @@ protected:
 };
 
 template<typename OutputIterator>
-void ColorScale::sample(int samples, OutputIterator it) const
+void ColorScale::sample(std::size_t samples, OutputIterator it) const
 {
     if (samples < 1) {
         return;
diff --git a/continuouscolorscale.cpp b/continuouscolorscale.cpp
index c729c6b..b7f2c93 100644
--- a/continuouscolorscale.cpp
+++ b/continuouscolorscale.cpp
@@ -9,8 +9,12 @@ ContinuousColorScale::ContinuousColorScale(std::initializer_list<QColor> colors)
 
 QColor ContinuousColorScale::color(float t) const
 {
-    if (t < m_min || t > m_max) {
-        return Color();
+    if (t <= m_min) {
+        return m_colors.front();
+    }
+
+    if (t >= m_max) {
+        return m_colors.back();
     }
 
     // normalize t
@@ -18,14 +22,10 @@ QColor ContinuousColorScale::color(float t) const
 
     // find which colors in the scale are adjacent to ours
     float step = 1.0 / m_colors.size();
-    int i = static_cast<int>(t / step);
-
-    if (i >= m_colors.size() - 1) {
-        return m_colors.back();
-    }
+    size_t i = static_cast<size_t>(t * m_colors.size());
 
     // normalize t between the two colors
-    int j = i + 1;
+    size_t j = i + 1;
     t = (t - i*step) / (j*step - i*step);
     return m_colors[i + round(t)];
 }
diff --git a/geometry.cpp b/geometry.cpp
index 101911e..1b514f0 100644
--- a/geometry.cpp
+++ b/geometry.cpp
@@ -4,13 +4,6 @@
 
 static const float PI = 3.1415f;
 
-Point2D::Point2D()
-{}
-
-void Point2D::set(float x, float y)
-{
-}
-
 RectF::RectF(float x, float y, float width, float height)
     : m_x(x)
     , m_y(y)
@@ -25,11 +18,11 @@ RectF::RectF(const RectF &other)
     , m_height(other.m_height)
 {}
 
-RectF::RectF(const Point2D &p1, const Point2D &p2)
-    : m_x(fmin(p1.x(), p2.x()))
-    , m_y(fmin(p1.y(), p2.y()))
-    , m_width(fabs(p1.x() - p2.x()))
-    , m_height(fabs(p1.y() - p2.y()))
+RectF::RectF(const vec2 &p1, const vec2 &p2)
+    : m_x(fmin(p1.x, p2.x))
+    , m_y(fmin(p1.y, p2.y))
+    , m_width(fabs(p1.x - p2.x))
+    , m_height(fabs(p1.y - p2.y))
 {}
 
 bool RectF::contains(float x, float y) const
@@ -38,6 +31,11 @@ bool RectF::contains(float x, float y) const
         && (y > m_y && y < m_y + m_height);
 }
 
+bool RectF::contains(const vec2 &p) const
+{
+    return contains(p.x, p.y);
+}
+
 bool RectF::intersects(const RectF &other) const
 {
     // TODO
@@ -52,7 +50,7 @@ int Geometry::vertexCount() const
     return 0;
 }
 
-Point2D *Geometry::vertexDataAsPoint2D()
+vec2 *Geometry::vertexDataAsPoint2D()
 {
     return nullptr;
 }
@@ -73,7 +71,7 @@ void updateCircleGeometry(Geometry *geometry, float diameter, float cx, float cy
     float x = diameter / 2;
     float y = 0;
 
-    Point2D *vertexData = geometry->vertexDataAsPoint2D();
+    vec2 *vertexData = geometry->vertexDataAsPoint2D();
     for (int i = 0; i < vertexCount; i++) {
         vertexData[i].set(x + cx, y + cy);
 
@@ -85,7 +83,7 @@ void updateCircleGeometry(Geometry *geometry, float diameter, float cx, float cy
 
 void updateRectGeometry(Geometry *geometry, float x, float y, float w, float h)
 {
-    Point2D *vertexData = geometry->vertexDataAsPoint2D();
+    vec2 *vertexData = geometry->vertexDataAsPoint2D();
 
     vertexData[0].set(x, y);
     vertexData[1].set(x + w, y);
@@ -99,7 +97,7 @@ void updateCrossHairGeometry(Geometry *geometry,
                              float thickness,
                              float length)
 {
-    Point2D *vertexData = geometry->vertexDataAsPoint2D();
+    vec2 *vertexData = geometry->vertexDataAsPoint2D();
     if (geometry->vertexCount() != 12) {
         return;
     }
diff --git a/geometry.h b/geometry.h
index 200782c..d7a10b8 100644
--- a/geometry.h
+++ b/geometry.h
@@ -1,17 +1,14 @@
 #ifndef GEOMETRY_H
 #define GEOMETRY_H
 
-class Point2D
+struct vec2
 {
-public:
-    Point2D();
+    float x, y;
 
-    void set(float, float);
-    float x() const { return m_x; }
-    float y() const { return m_y; }
-
-private:
-    float m_x, m_y;
+    void set(float _x, float _y) {
+        x = _x;
+        y = _y;
+    }
 };
 
 class RectF
@@ -19,7 +16,7 @@ class RectF
 public:
     RectF(float, float, float, float);
     RectF(const RectF &);
-    RectF(const Point2D &, const Point2D &);
+    RectF(const vec2 &, const vec2 &);
 
     float x() const { return m_x; }
     float y() const { return m_y; }
@@ -27,6 +24,7 @@ public:
     float height() const { return m_height; }
 
     bool contains(float x, float y) const;
+    bool contains(const vec2 &) const;
     bool intersects(const RectF &) const;
 
 private:
@@ -39,7 +37,7 @@ public:
     Geometry();
 
     int vertexCount() const;
-    Point2D *vertexDataAsPoint2D();
+    vec2 *vertexDataAsPoint2D();
 
 private:
 };
diff --git a/main.cpp b/main.cpp
index 5f76a74..17c6597 100644
--- a/main.cpp
+++ b/main.cpp
@@ -8,12 +8,34 @@
 
 #define GLAD_GL_IMPLEMENTATION
 #include <glad/gl.h>
+#undef GLAD_GL_IMPLEMENTATION
 #define GLFW_INCLUDE_NONE
 #include <GLFW/glfw3.h>
 
+#define NK_INCLUDE_FIXED_TYPES
+#define NK_INCLUDE_STANDARD_IO
+#define NK_INCLUDE_DEFAULT_ALLOCATOR
+#define NK_INCLUDE_VERTEX_BUFFER_OUTPUT
+#define NK_INCLUDE_FONT_BAKING
+#define NK_INCLUDE_DEFAULT_FONT
+#define NK_SQRT sqrt
+#define NK_SIN sinf
+#define NK_COS cosf
+// #define NK_STRTOD strtod
+// #define NK_DTOA dtoa
 #define NK_IMPLEMENTATION
 #include "nuklear.h"
 
+// device (NK) implementation details
+static const int MAX_VERTEX_MEMORY = 512 * 1024;
+static const int MAX_ELEMENT_MEMORY = 128 * 1024;
+
+#ifdef __APPLE__
+  #define NK_SHADER_VERSION "#version 150\n"
+#else
+  #define NK_SHADER_VERSION "#version 300 es\n"
+#endif
+
 #define STB_IMAGE_IMPLEMENTATION
 #include "stb_image.h"
 
@@ -36,6 +58,275 @@ static const int RNG_SEED = 123;
 static const int WINDOW_WIDTH = 1200;
 static const int WINDOW_HEIGHT = 800;
 
+struct media {
+    struct nk_font *font_14;
+    struct nk_font *font_18;
+    struct nk_font *font_20;
+    struct nk_font *font_22;
+
+    struct nk_image unchecked;
+    struct nk_image checked;
+    struct nk_image rocket;
+    struct nk_image cloud;
+    struct nk_image pen;
+    struct nk_image play;
+    struct nk_image pause;
+    struct nk_image stop;
+    struct nk_image prev;
+    struct nk_image next;
+    struct nk_image tools;
+    struct nk_image dir;
+    struct nk_image copy;
+    struct nk_image convert;
+    struct nk_image del;
+    struct nk_image edit;
+    struct nk_image images[9];
+    struct nk_image menu[6];
+};
+struct nk_glfw_vertex {
+    float position[2];
+    float uv[2];
+    nk_byte col[4];
+};
+
+struct device {
+    struct nk_buffer cmds;
+    struct nk_draw_null_texture tex_null;
+    GLuint vbo, vao, ebo;
+    GLuint prog;
+    GLuint vert_shdr;
+    GLuint frag_shdr;
+    GLint attrib_pos;
+    GLint attrib_uv;
+    GLint attrib_col;
+    GLint uniform_tex;
+    GLint uniform_proj;
+    GLuint font_tex;
+};
+
+static void
+ui_header(struct nk_context *ctx, struct media *media, const char *title)
+{
+    nk_style_set_font(ctx, &media->font_18->handle);
+    nk_layout_row_dynamic(ctx, 20, 1);
+    nk_label(ctx, title, NK_TEXT_LEFT);
+}
+
+static void
+ui_widget(struct nk_context *ctx, struct media *media, float height)
+{
+    static const float ratio[] = {0.15f, 0.85f};
+    nk_style_set_font(ctx, &media->font_22->handle);
+    nk_layout_row(ctx, NK_DYNAMIC, height, 2, ratio);
+    nk_spacing(ctx, 1);
+}
+
+static void
+device_init(struct device *dev)
+{
+    GLint status;
+    static const GLchar *vertex_shader =
+        NK_SHADER_VERSION
+        "uniform mat4 ProjMtx;\n"
+        "in vec2 Position;\n"
+        "in vec2 TexCoord;\n"
+        "in vec4 Color;\n"
+        "out vec2 Frag_UV;\n"
+        "out vec4 Frag_Color;\n"
+        "void main() {\n"
+        "   Frag_UV = TexCoord;\n"
+        "   Frag_Color = Color;\n"
+        "   gl_Position = ProjMtx * vec4(Position.xy, 0, 1);\n"
+        "}\n";
+    static const GLchar *fragment_shader =
+        NK_SHADER_VERSION
+        "precision mediump float;\n"
+        "uniform sampler2D Texture;\n"
+        "in vec2 Frag_UV;\n"
+        "in vec4 Frag_Color;\n"
+        "out vec4 Out_Color;\n"
+        "void main(){\n"
+        "   Out_Color = Frag_Color * texture(Texture, Frag_UV.st);\n"
+        "}\n";
+
+    nk_buffer_init_default(&dev->cmds);
+    dev->prog = glCreateProgram();
+    dev->vert_shdr = glCreateShader(GL_VERTEX_SHADER);
+    dev->frag_shdr = glCreateShader(GL_FRAGMENT_SHADER);
+    glShaderSource(dev->vert_shdr, 1, &vertex_shader, 0);
+    glShaderSource(dev->frag_shdr, 1, &fragment_shader, 0);
+    glCompileShader(dev->vert_shdr);
+    glCompileShader(dev->frag_shdr);
+    glGetShaderiv(dev->vert_shdr, GL_COMPILE_STATUS, &status);
+    assert(status == GL_TRUE);
+    glGetShaderiv(dev->frag_shdr, GL_COMPILE_STATUS, &status);
+    assert(status == GL_TRUE);
+    glAttachShader(dev->prog, dev->vert_shdr);
+    glAttachShader(dev->prog, dev->frag_shdr);
+    glLinkProgram(dev->prog);
+    glGetProgramiv(dev->prog, GL_LINK_STATUS, &status);
+    assert(status == GL_TRUE);
+
+    dev->uniform_tex = glGetUniformLocation(dev->prog, "Texture");
+    dev->uniform_proj = glGetUniformLocation(dev->prog, "ProjMtx");
+    dev->attrib_pos = glGetAttribLocation(dev->prog, "Position");
+    dev->attrib_uv = glGetAttribLocation(dev->prog, "TexCoord");
+    dev->attrib_col = glGetAttribLocation(dev->prog, "Color");
+
+    {
+        /* buffer setup */
+        GLsizei vs = sizeof(struct nk_glfw_vertex);
+        size_t vp = offsetof(struct nk_glfw_vertex, position);
+        size_t vt = offsetof(struct nk_glfw_vertex, uv);
+        size_t vc = offsetof(struct nk_glfw_vertex, col);
+
+        glGenBuffers(1, &dev->vbo);
+        glGenBuffers(1, &dev->ebo);
+        glGenVertexArrays(1, &dev->vao);
+
+        glBindVertexArray(dev->vao);
+        glBindBuffer(GL_ARRAY_BUFFER, dev->vbo);
+        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, dev->ebo);
+
+        glEnableVertexAttribArray((GLuint)dev->attrib_pos);
+        glEnableVertexAttribArray((GLuint)dev->attrib_uv);
+        glEnableVertexAttribArray((GLuint)dev->attrib_col);
+
+        glVertexAttribPointer((GLuint)dev->attrib_pos, 2, GL_FLOAT, GL_FALSE, vs, (void*)vp);
+        glVertexAttribPointer((GLuint)dev->attrib_uv, 2, GL_FLOAT, GL_FALSE, vs, (void*)vt);
+        glVertexAttribPointer((GLuint)dev->attrib_col, 4, GL_UNSIGNED_BYTE, GL_TRUE, vs, (void*)vc);
+    }
+
+    glBindTexture(GL_TEXTURE_2D, 0);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+    glBindVertexArray(0);
+}
+
+static void
+device_upload_atlas(struct device *dev, const void *image, int width, int height)
+{
+    glGenTextures(1, &dev->font_tex);
+    glBindTexture(GL_TEXTURE_2D, dev->font_tex);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, (GLsizei)width, (GLsizei)height, 0,
+                GL_RGBA, GL_UNSIGNED_BYTE, image);
+}
+
+static void
+device_shutdown(struct device *dev)
+{
+    glDetachShader(dev->prog, dev->vert_shdr);
+    glDetachShader(dev->prog, dev->frag_shdr);
+    glDeleteShader(dev->vert_shdr);
+    glDeleteShader(dev->frag_shdr);
+    glDeleteProgram(dev->prog);
+    glDeleteTextures(1, &dev->font_tex);
+    glDeleteBuffers(1, &dev->vbo);
+    glDeleteBuffers(1, &dev->ebo);
+    nk_buffer_free(&dev->cmds);
+}
+
+static void
+device_draw(struct device *dev, struct nk_context *ctx, int width, int height,
+    struct nk_vec2 scale, enum nk_anti_aliasing AA)
+{
+    GLfloat ortho[4][4] = {
+        {2.0f, 0.0f, 0.0f, 0.0f},
+        {0.0f,-2.0f, 0.0f, 0.0f},
+        {0.0f, 0.0f,-1.0f, 0.0f},
+        {-1.0f,1.0f, 0.0f, 1.0f},
+    };
+    ortho[0][0] /= (GLfloat)width;
+    ortho[1][1] /= (GLfloat)height;
+
+    /* setup global state */
+    glEnable(GL_BLEND);
+    glBlendEquation(GL_FUNC_ADD);
+    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+    glDisable(GL_CULL_FACE);
+    glDisable(GL_DEPTH_TEST);
+    glEnable(GL_SCISSOR_TEST);
+    glActiveTexture(GL_TEXTURE0);
+
+    /* setup program */
+    glUseProgram(dev->prog);
+    glUniform1i(dev->uniform_tex, 0);
+    glUniformMatrix4fv(dev->uniform_proj, 1, GL_FALSE, &ortho[0][0]);
+    {
+        /* convert from command queue into draw list and draw to screen */
+        const struct nk_draw_command *cmd;
+        void *vertices, *elements;
+        const nk_draw_index *offset = NULL;
+
+        /* allocate vertex and element buffer */
+        glBindVertexArray(dev->vao);
+        glBindBuffer(GL_ARRAY_BUFFER, dev->vbo);
+        glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, dev->ebo);
+
+        glBufferData(GL_ARRAY_BUFFER, MAX_VERTEX_MEMORY, NULL, GL_STREAM_DRAW);
+        glBufferData(GL_ELEMENT_ARRAY_BUFFER, MAX_ELEMENT_MEMORY, NULL, GL_STREAM_DRAW);
+
+        /* load draw vertices & elements directly into vertex + element buffer */
+        vertices = glMapBuffer(GL_ARRAY_BUFFER, GL_WRITE_ONLY);
+        elements = glMapBuffer(GL_ELEMENT_ARRAY_BUFFER, GL_WRITE_ONLY);
+        {
+            /* fill convert configuration */
+            struct nk_convert_config config;
+            static const struct nk_draw_vertex_layout_element vertex_layout[] = {
+                {NK_VERTEX_POSITION, NK_FORMAT_FLOAT, NK_OFFSETOF(struct nk_glfw_vertex, position)},
+                {NK_VERTEX_TEXCOORD, NK_FORMAT_FLOAT, NK_OFFSETOF(struct nk_glfw_vertex, uv)},
+                {NK_VERTEX_COLOR, NK_FORMAT_R8G8B8A8, NK_OFFSETOF(struct nk_glfw_vertex, col)},
+                {NK_VERTEX_LAYOUT_END}
+            };
+            NK_MEMSET(&config, 0, sizeof(config));
+            config.vertex_layout = vertex_layout;
+            config.vertex_size = sizeof(struct nk_glfw_vertex);
+            config.vertex_alignment = NK_ALIGNOF(struct nk_glfw_vertex);
+            config.tex_null = dev->tex_null;
+            config.circle_segment_count = 22;
+            config.curve_segment_count = 22;
+            config.arc_segment_count = 22;
+            config.global_alpha = 1.0f;
+            config.shape_AA = AA;
+            config.line_AA = AA;
+
+            /* setup buffers to load vertices and elements */
+            {struct nk_buffer vbuf, ebuf;
+            nk_buffer_init_fixed(&vbuf, vertices, MAX_VERTEX_MEMORY);
+            nk_buffer_init_fixed(&ebuf, elements, MAX_ELEMENT_MEMORY);
+            nk_convert(ctx, &dev->cmds, &vbuf, &ebuf, &config);}
+        }
+        glUnmapBuffer(GL_ARRAY_BUFFER);
+        glUnmapBuffer(GL_ELEMENT_ARRAY_BUFFER);
+
+        /* iterate over and execute each draw command */
+        nk_draw_foreach(cmd, ctx, &dev->cmds)
+        {
+            if (!cmd->elem_count) continue;
+            glBindTexture(GL_TEXTURE_2D, (GLuint)cmd->texture.id);
+            glScissor(
+                (GLint)(cmd->clip_rect.x * scale.x),
+                (GLint)((height - (GLint)(cmd->clip_rect.y + cmd->clip_rect.h)) * scale.y),
+                (GLint)(cmd->clip_rect.w * scale.x),
+                (GLint)(cmd->clip_rect.h * scale.y));
+            glDrawElements(GL_TRIANGLES, (GLsizei)cmd->elem_count, GL_UNSIGNED_SHORT, offset);
+            offset += cmd->elem_count;
+        }
+        nk_clear(ctx);
+        nk_buffer_clear(&dev->cmds);
+    }
+
+    /* default OpenGL state */
+    glUseProgram(0);
+    glBindBuffer(GL_ARRAY_BUFFER, 0);
+    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
+    glBindVertexArray(0);
+    glDisable(GL_BLEND);
+    glDisable(GL_SCISSOR_TEST);
+}
+
 class Main
 {
 public:
@@ -90,12 +381,27 @@ public:
     }
 
     enum ColorScaleType {
-        ColorScaleCategorical,
+        ColorScaleCategorical = 0,
         ColorScaleContinuous,
         ColorScaleDivergent,
         ColorScaleRainbow
     };
 
+    static constexpr const char *colormapItemNames[] = {
+        "Categorical",
+        "Continuous",
+        "Divergent",
+        "Rainbow",
+        nullptr,
+    };
+
+    static constexpr const char *metricItemNames[] = {
+        "Aggregate error",
+        "Ctrl. Pt. influence",
+        "Stress",
+        nullptr,
+    };
+
     void setCPColorScale(ColorScaleType colorScaleType) {
         float min = 0.0f;
         float max = 1.0f;
@@ -106,10 +412,10 @@ public:
         colorScaleCPs = getColorScale(colorScaleType);
         colorScaleCPs->setExtents(min, max);
 
-        cpPlot->setColorScale(colorScaleCPs);
-        cpBarChart->setColorScale(colorScaleCPs);
+        // cpBarChart->setColorScale(colorScaleCPs);
         cpColormap->setColorScale(colorScaleCPs);
         // bundlePlot->setColorScale(colorScaleCPs);
+        cpPlot->update();
     }
 
     void setRPColorScale(ColorScaleType colorScaleType) {
@@ -123,10 +429,10 @@ public:
         colorScaleRPs = getColorScale(colorScaleType);
         colorScaleRPs->setExtents(min, max);
 
-        rpPlot->setColorScale(colorScaleRPs);
-        splat->setColorScale(colorScaleRPs);
-        rpBarChart->setColorScale(colorScaleRPs);
+        // rpBarChart->setColorScale(colorScaleRPs);
         rpColormap->setColorScale(colorScaleRPs);
+        rpPlot->update();
+        splat->update();
     }
 
     // Pointers to visual components whose values are set in the main() function
@@ -238,6 +544,31 @@ private:
     std::string m_indicesSavePath, m_cpSavePath;
 };
 
+static void checkGLError(const char *msg) {
+    GLenum error = glGetError();
+    std::cout << msg << " ";
+    switch (error) {
+    case GL_NO_ERROR:
+        std::cout << "GL_NO_ERROR" << std::endl;
+        break;
+    case GL_INVALID_ENUM:
+        std::cout << "GL_INVALID_ENUM" << std::endl;
+        break;
+    case GL_INVALID_VALUE:
+        std::cout << "GL_INVALID_VALUE" << std::endl;
+        break;
+    case GL_INVALID_OPERATION:
+        std::cout << "GL_INVALID_OPERATION" << std::endl;
+        break;
+    case GL_INVALID_FRAMEBUFFER_OPERATION:
+        std::cout << "GL_INVALID_FRAMEBUFFER_OPERATION" << std::endl;
+        break;
+    case GL_OUT_OF_MEMORY:
+        std::cout << "GL_OUT_OF_MEMORY" << std::endl;
+        break;
+    }
+}
+
 arma::uvec extractCPs(const arma::mat &X)
 {
     int numCPs = (int) (3 * sqrt(X.n_rows));
@@ -247,7 +578,7 @@ arma::uvec extractCPs(const arma::mat &X)
     return indices.subvec(0, numCPs-1);
 }
 
-arma::mat standardize(const arma::mat &X)
+arma::mat standardise(const arma::mat &X)
 {
     arma::mat stdX = X;
     for (arma::uword j = 0; j < X.n_cols; j++) {
@@ -284,9 +615,9 @@ int main(int argc, char *argv[])
     int display_width = 0, display_height = 0;
 
     // GUI
-    // struct device device;
-    // struct nk_font_atlas atlas;
-    // struct media media;
+    struct device device;
+    struct nk_font_atlas atlas;
+    struct media media;
     struct nk_context ctx;
 
     NK_UNUSED(argc);
@@ -303,13 +634,13 @@ int main(int argc, char *argv[])
     glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
 #ifdef __APPLE__
     glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
-#endif
     glfwWindowHint(GLFW_SCALE_TO_MONITOR, GL_TRUE);
+#endif
     glfwWindowHint(GLFW_RED_BITS, 8);
     glfwWindowHint(GLFW_GREEN_BITS, 8);
     glfwWindowHint(GLFW_BLUE_BITS, 8);
     glfwWindowHint(GLFW_ALPHA_BITS, 8);
-    glfwWindowHint(GLFW_SAMPLES, 8);
+    glfwWindowHint(GLFW_SAMPLES, 4);
 
     win = glfwCreateWindow(WINDOW_WIDTH, WINDOW_HEIGHT, "pm", NULL, NULL);
     if (!win) {
@@ -342,15 +673,16 @@ int main(int argc, char *argv[])
         exit(1);
     }
 
+    glEnable(GL_MULTISAMPLE);
     glViewport(0, 0, WINDOW_WIDTH, WINDOW_HEIGHT);
 
     Main &m = Main::instance();
-    if (!m.loadDataset("wdbc-std.tbl")) {
+    if (argc != 2 || !m.loadDataset(argv[1])) {
         std::cerr << "Could not load data." << std::endl;
         return 1;
     }
 
-    const arma::mat &X = standardize(m.X());
+    const arma::mat &X = standardise(m.X());
 
     arma::arma_rng::set_seed(RNG_SEED);
     arma::uvec cpIndices;
@@ -371,7 +703,7 @@ int main(int argc, char *argv[])
 
     // TODO: maybe put this inside m->setCP() method; for now, will be outside
     // for more flexibility
-    Ys = arma::normalise(Ys);
+    // Ys = arma::normalise(Ys);
 
     // Sort indices (some operations become easier later)
     arma::uvec cpSortedIndices = arma::sort_index(cpIndices);
@@ -385,13 +717,23 @@ int main(int argc, char *argv[])
     ProjectionHistory history(X, cpIndices);
 
     // Visual components
+    Colormap cpColormap, rpColormap;
     Scatterplot cpPlot, rpPlot;
     VoronoiSplat splat;
-    Colormap cpColormap, rpColormap;
+    cpPlot.setSize(512, 512);
+    rpPlot.setSize(512, 512);
+    rpPlot.setGlyphSize(3.0f);
+    rpPlot.setColormap(rpColormap.texture());
+    splat.setSize(512, 512);
+    splat.setColormap(rpColormap.texture());
     BarChart cpBarChart, rpBarChart;
     m.cpPlot = &cpPlot;
     m.rpPlot = &rpPlot;
     m.splat = &splat;
+    m.cpColormap = &cpColormap;
+    m.rpColormap = &rpColormap;
+    m.cpBarChart = &cpBarChart;
+    m.rpBarChart = &rpBarChart;
 
     auto setCP = std::bind(&Main::setCP, &m, std::placeholders::_1);
     cpPlot.xyChanged.connect(setCP);
@@ -571,8 +913,8 @@ int main(int argc, char *argv[])
                                          std::max(m.colorScaleRPs->max(), (float) v.max()));
         }
 
-        m.splat->setColorScale(m.colorScaleRPs);
         m.rpColormap->setColorScale(m.colorScaleRPs);
+        // rpPlot.setColorData(v);
     });
     history.cpValuesChanged.connect([&](const arma::vec &v, bool rescale) {
         if (!m.colorScaleCPs || v.n_elem == 0) {
@@ -629,12 +971,40 @@ int main(int argc, char *argv[])
     // m->rpBarChart->setAcceptedMouseButtons(Qt::LeftButton | Qt::RightButton);
 
     m.setCPColorScale(Main::ColorScaleRainbow);
+    enum Main::ColorScaleType cpColorScaleType = Main::ColorScaleRainbow;
     m.setRPColorScale(Main::ColorScaleRainbow);
+    enum Main::ColorScaleType rpColorScaleType = Main::ColorScaleRainbow;
 
     // This sets the initial CP configuration, triggering all the necessary
     // signals to set up the helper objects and visual components
     manipulationHandler.setCP(Ys);
 
+
+    // GUI init
+    device_init(&device);
+    {
+        const void *image;
+        int w, h;
+        struct nk_font_config cfg = nk_font_config(0);
+        cfg.oversample_h = 3;
+        cfg.oversample_v = 2;
+        nk_font_atlas_init_default(&atlas);
+        nk_font_atlas_begin(&atlas);
+        media.font_14 = nk_font_atlas_add_from_file(&atlas, "assets/Roboto-Regular.ttf", 14.0f, &cfg);
+        media.font_18 = nk_font_atlas_add_from_file(&atlas, "assets/Roboto-Regular.ttf", 18.0f, &cfg);
+        media.font_20 = nk_font_atlas_add_from_file(&atlas, "assets/Roboto-Regular.ttf", 20.0f, &cfg);
+        media.font_22 = nk_font_atlas_add_from_file(&atlas, "assets/Roboto-Regular.ttf", 22.0f, &cfg);
+        image = nk_font_atlas_bake(&atlas, &w, &h, NK_FONT_ATLAS_RGBA32);
+        device_upload_atlas(&device, image, w, h);
+        nk_font_atlas_end(&atlas, nk_handle_id((int)device.font_tex), &device.tex_null);
+    }
+    nk_init_default(&ctx, &media.font_14->handle);
+
+    struct nk_image cpPlot_img = nk_image_id((int) cpPlot.texture());
+    struct nk_image rpPlot_img = nk_image_id((int) rpPlot.texture());
+    struct nk_image splat_img = nk_image_id((int) splat.texture());
+    struct nk_image cpColormap_img = nk_image_id((int) cpColormap.texture());
+    struct nk_image rpColormap_img = nk_image_id((int) rpColormap.texture());
     while (!glfwWindowShouldClose(win)) {
         struct nk_vec2 scale;
         glfwGetWindowSize(win, &width, &height);
@@ -673,19 +1043,210 @@ int main(int argc, char *argv[])
             int int_x = static_cast<int>(x);
             int int_y = static_cast<int>(y);
             nk_input_motion(&ctx, int_x, int_y);
-            nk_input_button(&ctx, NK_BUTTON_LEFT, int_x, int_y, glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_LEFT) == GLFW_PRESS);
-            nk_input_button(&ctx, NK_BUTTON_MIDDLE, int_x, int_y, glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_MIDDLE) == GLFW_PRESS);
-            nk_input_button(&ctx, NK_BUTTON_RIGHT, int_x, int_y, glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_RIGHT) == GLFW_PRESS);
+            nk_input_button(&ctx, NK_BUTTON_LEFT, int_x, int_y,
+                            glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_LEFT) == GLFW_PRESS);
+            nk_input_button(&ctx, NK_BUTTON_MIDDLE, int_x, int_y,
+                            glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_MIDDLE) == GLFW_PRESS);
+            nk_input_button(&ctx, NK_BUTTON_RIGHT, int_x, int_y,
+                            glfwGetMouseButton(win, GLFW_MOUSE_BUTTON_RIGHT) == GLFW_PRESS);
             nk_input_end(&ctx);
         }
 
+        nk_style_set_font(&ctx, &media.font_20->handle);
+        if (nk_begin(&ctx, "Control points", nk_rect(550, 20, 350, 300),
+                     NK_WINDOW_BORDER | NK_WINDOW_TITLE)) {
+            ui_header(&ctx, &media, "Glyphs");
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float glyphSize = cpPlot.glyphSize();
+            nk_property_float(&ctx, "Size:", 1.0f, &glyphSize, 20.0f, 1.0f, 0.1f);
+            if (glyphSize != cpPlot.glyphSize()) {
+                cpPlot.setGlyphSize(glyphSize);
+            }
+
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float opacity = 1.0f;
+            nk_property_float(&ctx, "Opacity:", 0.0f, &opacity, 1.0f, 0.1f, 0.01f);
+            // TODO
+
+            nk_layout_row_dynamic(&ctx, 20, 1);
+            nk_spacer(&ctx);
+
+            ui_header(&ctx, &media, "Colors");
+            nk_layout_row_begin(&ctx, NK_DYNAMIC, 30, 2);
+            nk_layout_row_push(&ctx, 0.3f);
+            nk_label(&ctx, "Scale:", NK_TEXT_RIGHT);
+            nk_layout_row_push(&ctx, 0.7f);
+            if (nk_combo_begin_label(&ctx, Main::colormapItemNames[cpColorScaleType],
+                                     nk_vec2(nk_widget_width(&ctx), 200))) {
+                nk_layout_row_dynamic(&ctx, 25, 1);
+                size_t i = 0;
+                for (const char * const *name = &Main::colormapItemNames[0];
+                     *name != nullptr; name++) {
+                    if (nk_combo_item_label(&ctx, *name, NK_TEXT_LEFT)) {
+                        enum Main::ColorScaleType type = static_cast<Main::ColorScaleType>(i);
+                        if (type != cpColorScaleType) {
+                            m.setCPColorScale(type);
+                            cpColorScaleType = type;
+                        }
+                    }
+                    i++;
+                }
+                nk_combo_end(&ctx);
+            }
+            nk_layout_row_end(&ctx);
+
+            nk_layout_row_begin(&ctx, NK_DYNAMIC, 30, 2);
+            nk_layout_row_push(&ctx, 0.3f);
+            nk_label(&ctx, "Map to:", NK_TEXT_RIGHT);
+            nk_layout_row_push(&ctx, 0.7f);
+            if (nk_combo_begin_label(&ctx, Main::metricItemNames[0],
+                                     nk_vec2(nk_widget_width(&ctx), 200))) {
+                nk_layout_row_dynamic(&ctx, 25, 1);
+                size_t i = 0;
+                for (const char * const *name = &Main::metricItemNames[0];
+                     *name != nullptr; name++) {
+                    if (nk_combo_item_label(&ctx, *name, NK_TEXT_LEFT)) {
+                        // TODO
+                    }
+                    i++;
+                }
+                nk_combo_end(&ctx);
+            }
+            nk_layout_row_end(&ctx);
+        }
+        nk_end(&ctx);
+
+        if (nk_begin(&ctx, "Regular points", nk_rect(550, 340, 350, 400),
+                     NK_WINDOW_BORDER | NK_WINDOW_TITLE)) {
+            ui_header(&ctx, &media, "Splat");
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float alpha = splat.alpha();
+            nk_property_float(&ctx, "Blur (alpha):", 5.0f, &alpha, 50.0f, 1.0f, 0.1f);
+            if (alpha != splat.alpha()) {
+                splat.setAlpha(alpha);
+            }
+
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float beta = splat.beta();
+            nk_property_float(&ctx, "Radius (beta):", 5.0f, &beta, 50.0f, 1.0f, 0.1f);
+            if (beta != splat.beta()) {
+                splat.setBeta(beta);
+            }
+
+            nk_layout_row_dynamic(&ctx, 20, 1);
+            nk_spacer(&ctx);
+
+            ui_header(&ctx, &media, "Glyphs");
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float glyphSize = rpPlot.glyphSize();
+            nk_property_float(&ctx, "Size:", 1.0f, &glyphSize, 20.0f, 1.0f, 0.1f);
+            if (glyphSize != rpPlot.glyphSize()) {
+                rpPlot.setGlyphSize(glyphSize);
+            }
+
+            nk_layout_row_dynamic(&ctx, 30, 1);
+            float opacity = 1.0f;
+            nk_property_float(&ctx, "Opacity:", 0.0f, &opacity, 1.0f, 0.1f, 0.01f);
+            // TODO
+
+            nk_layout_row_dynamic(&ctx, 20, 1);
+            nk_spacer(&ctx);
+
+            ui_header(&ctx, &media, "Colors");
+            nk_layout_row_begin(&ctx, NK_DYNAMIC, 30, 2);
+            nk_layout_row_push(&ctx, 0.3f);
+            nk_label(&ctx, "Scale:", NK_TEXT_RIGHT);
+            nk_layout_row_push(&ctx, 0.7f);
+            if (nk_combo_begin_label(&ctx, Main::colormapItemNames[rpColorScaleType],
+                                     nk_vec2(nk_widget_width(&ctx), 200))) {
+                nk_layout_row_dynamic(&ctx, 25, 1);
+                size_t i = 0;
+                for (const char * const *name = &Main::colormapItemNames[0];
+                     *name != nullptr; name++) {
+                    if (nk_combo_item_label(&ctx, *name, NK_TEXT_LEFT)) {
+                        enum Main::ColorScaleType type = static_cast<Main::ColorScaleType>(i);
+                        if (type != rpColorScaleType) {
+                            m.setRPColorScale(type);
+                            rpColorScaleType = type;
+                        }
+                    }
+                    i++;
+                }
+                nk_combo_end(&ctx);
+            }
+            nk_layout_row_end(&ctx);
+
+            nk_layout_row_begin(&ctx, NK_DYNAMIC, 30, 2);
+            nk_layout_row_push(&ctx, 0.3f);
+            nk_label(&ctx, "Map to:", NK_TEXT_RIGHT);
+            nk_layout_row_push(&ctx, 0.7f);
+            if (nk_combo_begin_label(&ctx, Main::metricItemNames[0],
+                                     nk_vec2(nk_widget_width(&ctx), 200))) {
+                nk_layout_row_dynamic(&ctx, 25, 1);
+                size_t i = 0;
+                for (const char * const *name = &Main::metricItemNames[0];
+                     *name != nullptr; name++) {
+                    if (nk_combo_item_label(&ctx, *name, NK_TEXT_LEFT)) {
+                        // TODO
+                    }
+                    i++;
+                }
+                nk_combo_end(&ctx);
+            }
+            nk_layout_row_end(&ctx);
+        }
+        nk_end(&ctx);
+
+        if (nk_begin(&ctx, "Scatterplot", nk_rect(20, 20, splat.width(), splat.height()),
+                     NK_WINDOW_BORDER | NK_WINDOW_NO_SCROLLBAR)) {
+            struct nk_command_buffer *canvas = nk_window_get_canvas(&ctx);
+            struct nk_rect region{0.0f, 0.0f,
+                                  static_cast<float>(splat.width()),
+                                  static_cast<float>(splat.height())};
+
+            // Render components to their textures
+            splat.draw();
+            rpPlot.draw();
+            cpPlot.draw();
+
+            // Add white background rect and draw textures on top
+            nk_layout_space_begin(&ctx, NK_STATIC, region.h, 4);
+            nk_layout_space_push(&ctx, region);
+            nk_fill_rect(canvas, region, 0.0f, nk_rgba(255, 255, 255, 255));
+            nk_layout_space_push(&ctx, region);
+            nk_image(&ctx, splat_img);
+            nk_layout_space_push(&ctx, region);
+            nk_image(&ctx, rpPlot_img);
+            nk_layout_space_push(&ctx, region);
+            nk_image(&ctx, cpPlot_img);
+            nk_layout_space_end(&ctx);
+        }
+        nk_end(&ctx);
+
+        if (nk_begin(&ctx, "Colormap", nk_rect(20, 40 + splat.height(), splat.width(), 120),
+                     NK_WINDOW_BORDER | NK_WINDOW_TITLE | NK_WINDOW_NO_SCROLLBAR)) {
+            // nk_widget(&bounds, &ctx);
+            ui_header(&ctx, &media, "Control points");
+            nk_layout_row_dynamic(&ctx, 10, 1);
+            nk_image(&ctx, cpColormap_img);
+
+            ui_header(&ctx, &media, "Regular points");
+            nk_layout_row_dynamic(&ctx, 10, 1);
+            nk_image(&ctx, rpColormap_img);
+        }
+        nk_end(&ctx);
+
         glViewport(0, 0, display_width, display_height);
-        glClear(GL_COLOR_BUFFER_BIT);
+        // glDisable(GL_BLEND);
         glClearColor(0.3f, 0.3f, 0.3f, 1.0f);
-        // device_draw(&device, &ctx, width, height, scale, NK_ANTI_ALIASING_ON);
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+        device_draw(&device, &ctx, width, height, scale, NK_ANTI_ALIASING_ON);
+
         glfwSwapBuffers(win);
+        // nk_clear(&ctx);
     }
 
+    nk_font_atlas_clear(&atlas);
     nk_free(&ctx);
     glfwTerminate();
     return 0;
diff --git a/manifest.scm b/manifest.scm
index 2493f02..83035ed 100644
--- a/manifest.scm
+++ b/manifest.scm
@@ -1,6 +1,7 @@
 (use-modules (guix packages)
              (gnu packages cmake)
              (gnu packages commencement)
+             (gnu packages cpp)
              (gnu packages gdb)
              (gnu packages gl)
              (gnu packages maths)
@@ -11,6 +12,7 @@
 (packages->manifest
  (list
   armadillo
+  ccls
   cmake
   egl-wayland
   gcc-toolchain
diff --git a/scale.cpp b/scale.cpp
new file mode 100644
index 0000000..ba0cff0
--- /dev/null
+++ b/scale.cpp
@@ -0,0 +1,38 @@
+#include "scale.h"
+
+#include <iostream>
+
+void updateTransform4x4(const LinearScale<float> &sx,
+                        const LinearScale<float> &sy,
+                        float transform[4][4])
+{
+    float fsx = 2.0f * sx.slope();
+    float ftx = 2.0f * sx.offset() - 1.0f;
+    float fsy = 2.0f * sy.slope();
+    float fty = 2.0f * sy.offset() - 1.0f;
+
+    // The transform matrix should be this, but we have it transposed
+    // in memory:
+    //
+    // [   sx  0.0f  0.0f    tx ]
+    // [ 0.0f    sy  0.0f    ty ]
+    // [ 0.0f  0.0f  0.0f  0.0f ]
+    // [ 0.0f  0.0f  0.0f  1.0f ]
+    //
+    // It is already initialized with 0s and bottom-right 1
+    transform[0][0] = fsx;
+    transform[1][1] = fsy;
+    transform[3][0] = ftx;
+    transform[3][1] = fty;
+}
+
+void updateTransform4x4(LinearScale<float> &sx,
+                        LinearScale<float> &sy,
+                        float offsetX,
+                        float offsetY,
+                        float transform[4][4])
+{
+    sx.setRange(offsetX, 1.0f - offsetX);
+    sy.setRange(1.0f - offsetY, offsetY);
+    updateTransform4x4(sx, sy, transform);
+}
diff --git a/scale.h b/scale.h
index 4c3f341..b383f19 100644
--- a/scale.h
+++ b/scale.h
@@ -103,4 +103,14 @@ private:
     T m_transformSlope;
 };
 
+void updateTransform4x4(const LinearScale<float> &sx,
+                        const LinearScale<float> &sy,
+                        float transform[4][4]);
+
+void updateTransform4x4(LinearScale<float> &sx,
+                        LinearScale<float> &sy,
+                        float offsetX,
+                        float offsetY,
+                        float transform[4][4]);
+
 #endif // SCALE_H
diff --git a/scatterplot.cpp b/scatterplot.cpp
index fc31daa..a7194e8 100644
--- a/scatterplot.cpp
+++ b/scatterplot.cpp
@@ -3,8 +3,8 @@
 #include <algorithm>
 #include <cmath>
 
+#include "colormap.h"
 #include "continuouscolorscale.h"
-#include "geometry.h"
 
 // Glyphs settings
 static const Color DEFAULT_GLYPH_COLOR(255, 255, 255);
@@ -31,28 +31,197 @@ static const Color SELECTION_COLOR(128, 128, 128, 96);
 // static const Qt::MouseButton NORMAL_BUTTON = Qt::LeftButton;
 // static const Qt::MouseButton SPECIAL_BUTTON = Qt::RightButton;
 
+static const char *shaderVertex = R"EOF(#version 330
+
+precision mediump float;
+
+uniform float colormask;
+uniform mat4 transform;
+uniform sampler2D colormap;
+
+layout (location = 0) in vec2 pos;
+layout (location = 1) in float value;
+// layout (location = 2) in float opacity;
+
+out VS_OUT {
+    vec4 color;
+} vs_out;
+
+vec3 getRGB(float value) {
+    return texture(colormap, vec2(mix(0.005, 0.995, value), 0)).rgb;
+}
+
+void main()
+{
+    vs_out.color = vec4(mix(vec3(1.0), getRGB(value), colormask), 1.0);
+    // vs_out.color = vec4(1.0);
+    gl_Position = transform * vec4(pos.xy, 0.0, 1.0);
+}
+)EOF";
+
+static const char *shaderFragment = R"EOF(#version 330
+in vec4 color;
+
+out vec4 FragColor;
+
+void main()
+{
+    FragColor = color;
+    // FragColor = vec4(1.0);
+}
+)EOF";
+
+static const char *shaderGeometry = R"EOF(#version 330
+layout (points) in;
+layout (triangle_strip, max_vertices = 16) out;
+
+uniform float size;
+
+in VS_OUT {
+    vec4 color;
+} gs_in[];
+
+out vec4 color;
+
+#define PI 3.1415926538
+
+void main() {
+    int num = 8;
+    float theta = 2 * PI / num;
+    float c = cos(theta);
+    float s = sin(theta);
+    vec4 p = vec4(size, 0.0, 0.0, 0.0);
+    // gs_in[0] as there is only one input vertex
+    color = gs_in[0].color;
+    for (int i = 0; i < num; i++) {
+        gl_Position = gl_in[0].gl_Position + p;
+        EmitVertex();
+
+        gl_Position = gl_in[0].gl_Position + vec4(p.x, -p.y, 0.0, 0.0);
+        EmitVertex();
+
+        p.xy = vec2(c*p.x - s*p.y, s*p.x + c*p.y);
+    }
+    EndPrimitive();
+}
+)EOF";
+
+static const char *shaderOutlineFragment = R"EOF(#version 330
+in vec4 color;
+
+out vec4 FragColor;
+
+void main()
+{
+    FragColor = vec4(0.0, 0.0, 0.0, 1.0);
+}
+)EOF";
+
+static const char *shaderOutlineGeometry = R"EOF(#version 330
+layout (points) in;
+layout (line_strip, max_vertices = 16) out;
+
+uniform float size;
+
+in VS_OUT {
+    vec4 color;
+} gs_in[];
+
+out vec4 color;
+
+#define PI 3.1415926538
+
+void main() {
+    int num = 15;
+    float theta = 2 * PI / num;
+    float c = cos(theta);
+    float s = sin(theta);
+    vec4 p = vec4(size, 0.0, 0.0, 0.0);
+    // gs_in[0] as there is only one input vertex
+    color = gs_in[0].color;
+    for (int i = 0; i < num; i++) {
+        gl_Position = gl_in[0].gl_Position + p;
+        EmitVertex();
+
+        p.xy = vec2(c*p.x - s*p.y, s*p.x + c*p.y);
+    }
+    gl_Position = gl_in[0].gl_Position + p;
+    EmitVertex();
+    EndPrimitive();
+}
+)EOF";
+
 Scatterplot::Scatterplot()
     : m_glyphSize(DEFAULT_GLYPH_SIZE)
     , m_colorScale(0)
-    , m_autoScale(true)
+    , m_autoScale(false)
     , m_sx(0, 1, 0, 1)
     , m_sy(0, 1, 0, 1)
     , m_anySelected(false)
     , m_brushedItem(-1)
     , m_interactionState(StateNone)
     , m_dragEnabled(false)
+    , m_redraw(false)
     , m_shouldUpdateGeometry(false)
     , m_shouldUpdateMaterials(false)
 {
+    glGenFramebuffers(1, &m_FBO);
+
+    glGenBuffers(1, &m_pointsVBO);
+    glGenBuffers(1, &m_valuesVBO);
+    glGenVertexArrays(1, &m_VAO);
+    glBindVertexArray(m_VAO);
+    glBindBuffer(GL_ARRAY_BUFFER, m_pointsVBO);
+    glEnableVertexAttribArray(0);
+    glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, nullptr);
+    // glVertexAttribDivisor(0, 1);
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_valuesVBO);
+    glEnableVertexAttribArray(1);
+    glVertexAttribPointer(1, 1, GL_FLOAT, GL_FALSE, 0, nullptr);
+    glVertexAttribDivisor(1, 1);
+    glBindVertexArray(0);
+
+    glGenTextures(1, &m_outTex);
+
+    m_shader = std::make_unique<Shader>(
+        shaderVertex,
+        shaderFragment,
+        shaderGeometry);
+
+    m_shaderOutline = std::make_unique<Shader>(
+        shaderVertex,
+        shaderOutlineFragment,
+        shaderOutlineGeometry);
+
+    std::fill(&m_transform[0][0], &m_transform[0][0] + 16, 0.0f);
+    m_transform[3][3] = 1.0f;
 }
 
-void Scatterplot::setColorScale(std::shared_ptr<const ColorScale> colorScale)
+Scatterplot::~Scatterplot()
 {
-    m_colorScale = colorScale;
-    if (m_colorData.n_elem > 0) {
-        m_shouldUpdateMaterials = true;
-        update();
-    }
+    glDeleteFramebuffers(1, &m_FBO);
+}
+
+void Scatterplot::setSize(size_t width, size_t height)
+{
+    m_width = width;
+    m_height = height;
+
+    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);
+
+    update();
+}
+
+void Scatterplot::setColormap(GLuint texture)
+{
+    m_colormapTex = texture;
+    update();
 }
 
 arma::mat Scatterplot::XY() const
@@ -122,8 +291,6 @@ void Scatterplot::setScale(const LinearScale<float> &sx, const LinearScale<float
     scaleChanged(m_sx, m_sy);
 
     updateQuadTree();
-
-    m_shouldUpdateGeometry = true;
     update();
 }
 
@@ -150,13 +317,115 @@ void Scatterplot::setGlyphSize(float glyphSize)
 
     m_glyphSize = glyphSize;
     glyphSizeChanged(m_glyphSize);
-
-    m_shouldUpdateGeometry = true;
     update();
 }
 
 void Scatterplot::update()
 {
+    m_redraw = true;
+
+    if (m_shouldUpdateGeometry) {
+        updateGeometry();
+    }
+    if (m_shouldUpdateMaterials) {
+        updateMaterials();
+    }
+    updateTransform();
+}
+
+void Scatterplot::updateGeometry()
+{
+    if (m_xy.n_rows < 1 || m_xy.n_cols != 2) {
+        return;
+    }
+
+    if (m_colorData.n_elem > 0 && m_colorData.n_elem != m_xy.n_rows) {
+        // Old values are no longer valid, clean up
+        m_colorData.fill(0.0f);
+    }
+
+    // Copy 'points' to internal data structure(s)
+    m_points.resize(m_xy.n_rows);
+    const double *col_x = m_xy.colptr(0);
+    const double *col_y = m_xy.colptr(1);
+    for (unsigned i = 0; i < m_xy.n_rows; i++) {
+        m_points[i].set(col_x[i], col_y[i]);
+    }
+
+    glBindBuffer(GL_ARRAY_BUFFER, m_pointsVBO);
+    glBufferData(GL_ARRAY_BUFFER, m_points.size() * sizeof(vec2),
+                 m_points.data(), GL_DYNAMIC_DRAW);
+
+    m_shouldUpdateGeometry = false;
+}
+
+void Scatterplot::updateMaterials()
+{
+    glBindBuffer(GL_ARRAY_BUFFER, m_valuesVBO);
+    glBufferData(GL_ARRAY_BUFFER, m_colorData.n_elem * sizeof(float),
+                 m_colorData.memptr(), GL_DYNAMIC_DRAW);
+
+    m_shouldUpdateMaterials = false;
+}
+
+void Scatterplot::updateTransform()
+{
+    float padding = Scatterplot::PADDING;
+    float offsetX = padding / m_width, offsetY = padding / m_height;
+    updateTransform4x4(m_sx, m_sy, offsetX, offsetY, m_transform);
+}
+
+void Scatterplot::draw()
+{
+    if (!m_redraw) {
+        return;
+    }
+
+    int originalFBO;
+    glGetIntegerv(GL_FRAMEBUFFER_BINDING, &originalFBO);
+
+    glBindFramebuffer(GL_FRAMEBUFFER, m_FBO);
+    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
+                           GL_TEXTURE_2D, m_outTex, 0);
+
+    glViewport(0, 0, m_width, m_height);
+    m_shader->use();
+    m_shader->setUniform("transform", m_transform);
+    m_shader->setUniform("size", m_glyphSize / m_width);
+    m_shader->setUniform("colormask", m_colorData.n_elem == 0 ? 0.0f : 1.0f);
+
+    glActiveTexture(GL_TEXTURE0);
+    glBindTexture(GL_TEXTURE_2D, m_colormapTex);
+    m_shader->setUniform("colormap", 0);
+
+    glEnable(GL_BLEND);
+    // glBlendFuncSeparate(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA, GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
+    glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
+    glClear(GL_COLOR_BUFFER_BIT);
+
+    glBindVertexArray(m_VAO);
+    glDrawArrays(GL_POINTS, 0, m_points.size());
+    m_shader->release();
+
+    // Swap out shaders to draw outline
+    m_shaderOutline->use();
+    m_shaderOutline->setUniform("transform", m_transform);
+    m_shaderOutline->setUniform("size", m_glyphSize / m_width);
+    m_shaderOutline->setUniform("colormap", 0);
+    glDrawArrays(GL_POINTS, 0, m_points.size());
+
+    glBindVertexArray(0);
+    m_shaderOutline->release();
+
+    /*
+    GLenum status = glCheckFramebufferStatus(GL_FRAMEBUFFER);
+    if (status != GL_FRAMEBUFFER_COMPLETE) {
+        checkGLError("draw():");
+    }
+    */
+
+    glBindFramebuffer(GL_FRAMEBUFFER, originalFBO);
+    m_redraw = false;
 }
 
 /*
@@ -590,8 +859,8 @@ void Scatterplot::applyManipulation()
     m_sy.inverse();
     m_sx.inverse();
 
-    float tx = m_dragCurrentPos.x() - m_dragOriginPos.x();
-    float ty = m_dragCurrentPos.y() - m_dragOriginPos.y();
+    float tx = m_dragCurrentPos.x - m_dragOriginPos.x;
+    float ty = m_dragCurrentPos.y - m_dragOriginPos.y;
 
     for (std::vector<bool>::size_type i = 0; i < m_selection.size(); i++) {
         if (m_selection[i]) {
@@ -609,10 +878,10 @@ void Scatterplot::applyManipulation()
 
 void Scatterplot::updateQuadTree()
 {
-    m_sx.setRange(PADDING, width() - PADDING);
-    m_sy.setRange(height() - PADDING, PADDING);
+    m_sx.setRange(PADDING, m_width - PADDING);
+    m_sy.setRange(m_height - PADDING, PADDING);
 
-    m_quadtree.reset(new QuadTree(RectF(x(), y(), width(), height())));
+    m_quadtree.reset(new QuadTree(RectF(0.0f, 0.0f, m_width, m_height)));
     for (arma::uword i = 0; i < m_xy.n_rows; i++) {
         const arma::rowvec &row = m_xy.row(i);
         m_quadtree->insert(m_sx(row[0]), m_sy(row[1]), (int) i);
diff --git a/scatterplot.h b/scatterplot.h
index 8286b7d..5be1647 100644
--- a/scatterplot.h
+++ b/scatterplot.h
@@ -5,12 +5,14 @@
 #include <vector>
 
 #include <armadillo>
+#include <glad/gl.h>
 #include <nod.hpp>
 
 #include "colorscale.h"
 #include "quadtree.h"
 #include "geometry.h"
 #include "scale.h"
+#include "shader.h"
 
 class Scatterplot
 {
@@ -18,15 +20,15 @@ public:
     static const int PADDING = 20;
 
     Scatterplot();
+    ~Scatterplot();
 
     arma::mat XY() const;
-    void setColorScale(std::shared_ptr<const ColorScale> colorScale);
+    void setColormap(GLuint texture);
     void setAutoScale(bool autoScale);
 
-    float x() const { return m_x; }
-    float y() const { return m_y; }
-    float width() const { return m_width; }
-    float height() const { return m_height; }
+    size_t width() const { return m_width; }
+    size_t height() const { return m_height; }
+    void setSize(size_t, size_t);
 
     float glyphSize() const { return m_glyphSize; }
     void setGlyphSize(float glyphSize);
@@ -34,6 +36,8 @@ public:
     void setDragEnabled(bool enabled) { m_dragEnabled = enabled; }
     bool isDragEnabled() const { return m_dragEnabled; }
 
+    GLuint texture() const { return m_outTex; }
+
     nod::signal<void(const arma::mat &)> xyChanged, xyInteractivelyChanged;
     nod::signal<void(const arma::vec &)> colorDataChanged, opacityDataChanged;
     nod::signal<void(const std::vector<bool> &)> selectionChanged, selectionInteractivelyChanged;
@@ -48,6 +52,7 @@ public:
     void setSelection(const std::vector<bool> &selection);
     void brushItem(int item);
     void update();
+    void draw();
 
 protected:
     // QSGNode *updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *);
@@ -64,6 +69,9 @@ private:
     // QSGNode *newGlyphTree();
 
     void applyManipulation();
+    void updateGeometry();
+    void updateMaterials();
+    void updateTransform();
     // void updateGlyphs(QSGNode *node);
     // void updateBrush(QSGNode *node);
 
@@ -73,7 +81,7 @@ private:
     arma::vec m_opacityData;
 
     // Visuals
-    float m_x, m_y, m_width, m_height;
+    size_t m_width, m_height;
     float m_glyphSize;
     std::shared_ptr<const ColorScale> m_colorScale;
 
@@ -96,9 +104,14 @@ private:
     } m_interactionState;
     bool m_dragEnabled;
 
-    Point2D m_dragOriginPos, m_dragCurrentPos;
+    std::unique_ptr<Shader> m_shader, m_shaderOutline;
+    GLuint m_VAO, m_pointsVBO, m_valuesVBO;
+    GLuint m_FBO, m_colormapTex, m_outTex;
+    GLfloat m_transform[4][4];
+    std::vector<vec2> m_points;
 
-    bool m_shouldUpdateGeometry, m_shouldUpdateMaterials;
+    vec2 m_dragOriginPos, m_dragCurrentPos;
+    bool m_redraw, m_shouldUpdateGeometry, m_shouldUpdateMaterials;
 
     std::unique_ptr<QuadTree> m_quadtree;
     void updateQuadTree();
diff --git a/shader.cpp b/shader.cpp
new file mode 100644
index 0000000..0a1268e
--- /dev/null
+++ b/shader.cpp
@@ -0,0 +1,137 @@
+#include "shader.h"
+
+#include <iostream>
+
+static int setupShader(GLenum shaderType, const char *src)
+{
+    GLuint shader = glCreateShader(shaderType);
+    glShaderSource(shader, 1, &src, nullptr);
+    glCompileShader(shader);
+
+    int status;
+    glGetShaderiv(shader, GL_COMPILE_STATUS, &status);
+    if (status == GL_FALSE) {
+        char buf[1024];
+        glGetShaderInfoLog(shader, sizeof buf, nullptr, buf);
+        std::cerr << "Error: "
+                  << (shaderType == GL_VERTEX_SHADER ? "vertex" : "fragment")
+                  << " shader: compilation failed: "
+                  << buf
+                  << std::endl;
+    }
+    return shader;
+}
+
+Shader::Shader(const char *vertexSrc, const char *fragmentSrc, const char *geometrySrc)
+{
+    GLuint vertex = setupShader(GL_VERTEX_SHADER, vertexSrc);
+    GLuint fragment = setupShader(GL_FRAGMENT_SHADER, fragmentSrc);
+    GLuint geometry;
+    if (geometrySrc != nullptr) {
+        geometry = setupShader(GL_GEOMETRY_SHADER, geometrySrc);
+    }
+
+    int status;
+    m_program = glCreateProgram();
+    glAttachShader(m_program, vertex);
+    glAttachShader(m_program, fragment);
+    if (geometrySrc != nullptr) {
+        glAttachShader(m_program, geometry);
+    }
+    glLinkProgram(m_program);
+    glGetProgramiv(m_program, GL_LINK_STATUS, &status);
+    if (status == GL_FALSE) {
+        char buf[1024];
+        glGetProgramInfoLog(m_program, sizeof buf, nullptr, buf);
+        std::cerr << "Error: shader linking failed: "
+                  << buf
+                  << std::endl;
+    }
+
+    if (geometrySrc != nullptr) {
+        glDeleteShader(geometry);
+    }
+    glDeleteShader(fragment);
+    glDeleteShader(vertex);
+}
+
+Shader::~Shader()
+{
+    glDeleteProgram(m_program);
+}
+
+void Shader::use() const
+{
+    glUseProgram(m_program);
+}
+
+void Shader::release() const
+{
+    glUseProgram(0);
+}
+
+GLint Shader::uniformLocation(const char *name) const
+{
+    return glGetUniformLocation(m_program, name);
+}
+
+void Shader::setUniform(GLint location, GLint value) const
+{
+    glUniform1i(location, value);
+}
+
+void Shader::setUniform(GLint location, GLfloat value) const
+{
+    glUniform1f(location, value);
+}
+
+void Shader::setUniform(GLint location, GLfloat mat[4][4]) const
+{
+    // GL_FALSE for column-major
+    glUniformMatrix4fv(location, 1, GL_FALSE, &mat[0][0]);
+}
+
+void Shader::setUniform1dArray(GLint location, const GLfloat *values, int count) const
+{
+    glUniform1fv(location, count, values);
+}
+
+void Shader::setUniform2dArray(GLint location, const GLfloat *values, int count) const
+{
+    glUniform2fv(location, count, values);
+}
+
+void Shader::setUniform3dArray(GLint location, const GLfloat *values, int count) const
+{
+    glUniform3fv(location, count, values);
+}
+
+void Shader::setUniform(const char *name, GLint value) const
+{
+    setUniform(uniformLocation(name), value);
+}
+
+void Shader::setUniform(const char *name, GLfloat value) const
+{
+    setUniform(uniformLocation(name), value);
+}
+
+void Shader::setUniform1dArray(const char *name, const GLfloat *values, int count) const
+{
+    setUniform1dArray(uniformLocation(name), values, count);
+}
+
+void Shader::setUniform2dArray(const char *name, const GLfloat *values, int count) const
+{
+    setUniform2dArray(uniformLocation(name), values, count);
+}
+
+void Shader::setUniform3dArray(const char *name, const GLfloat *values, int count) const
+{
+    setUniform3dArray(uniformLocation(name), values, count);
+}
+
+void Shader::setUniform(const char *name, GLfloat mat[4][4]) const
+{
+    setUniform(uniformLocation(name), mat);
+}
diff --git a/shader.h b/shader.h
new file mode 100644
index 0000000..feb5952
--- /dev/null
+++ b/shader.h
@@ -0,0 +1,37 @@
+#ifndef SHADER_H
+#define SHADER_H
+
+#include <glad/gl.h>
+
+class Shader
+{
+public:
+    Shader(const char *vertexSrc, const char *fragmentSrc, const char *geometrySrc = nullptr);
+    ~Shader();
+
+    GLuint program() const { return m_program; }
+
+    GLint uniformLocation(const char *name) const;
+
+    void setUniform(GLint location, GLint value) const;
+    void setUniform(GLint location, GLfloat value) const;
+    void setUniform(GLint location, GLfloat mat[4][4]) const;
+    void setUniform1dArray(GLint location, const GLfloat *values, int count) const;
+    void setUniform2dArray(GLint location, const GLfloat *values, int count) const;
+    void setUniform3dArray(GLint location, const GLfloat *values, int count) const;
+
+    void setUniform(const char *name, GLint value) const;
+    void setUniform(const char *name, GLfloat value) const;
+    void setUniform(const char *name, GLfloat mat[4][4]) const;
+    void setUniform1dArray(const char *name, const GLfloat *values, int count) const;
+    void setUniform2dArray(const char *name, const GLfloat *values, int count) const;
+    void setUniform3dArray(const char *name, const GLfloat *values, int count) const;
+
+    void use() const;
+    void release() const;
+
+private:
+    GLuint m_program;
+};
+
+#endif // SHADER_H
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
diff --git a/voronoisplat.h b/voronoisplat.h
index 894b6ca..bab3fc0 100644
--- a/voronoisplat.h
+++ b/voronoisplat.h
@@ -4,20 +4,25 @@
 #include <memory>
 
 #include <armadillo>
+#include <glad/gl.h>
 #include <nod.hpp>
 
 #include "colorscale.h"
+#include "geometry.h"
 #include "scale.h"
+#include "shader.h"
 
 class VoronoiSplat
 {
 public:
     VoronoiSplat();
+    ~VoronoiSplat();
 
     // Renderer *createRenderer() const;
     void update();
+    void draw();
 
-    const std::vector<float> &sites() const      { return m_sites; }
+    const std::vector<vec2> &sites() const      { return m_sites; }
     const std::vector<float> &values() const     { return m_values; }
     const std::vector<float> &colorScale() const { return m_cmap; }
     LinearScale<float> scaleX() const { return m_sx; }
@@ -25,6 +30,11 @@ public:
     float alpha() const { return m_alpha; }
     float beta() const  { return m_beta; }
 
+    size_t width() const { return m_width; }
+    size_t height() const { return m_height; }
+
+    GLuint texture() { return m_outTex; }
+
     void setSitesChanged(bool sitesChanged) {
         m_sitesChanged = sitesChanged;
     }
@@ -37,18 +47,20 @@ public:
 
     nod::signal<void(const arma::mat &)> sitesChanged;
     nod::signal<void(const arma::vec &)> valuesChanged;
-    nod::signal<void(std::shared_ptr<const ColorScale>)> colorScaleChanged;
+    nod::signal<void(GLuint)> colormapChanged;
     nod::signal<void(const LinearScale<float> &, const LinearScale<float> &)> scaleChanged;
     nod::signal<void(float)> alphaChanged, betaChanged;
 
+    void setSize(size_t width, size_t height);
+
     // 'points' should be a 2D points matrix (each point in a row)
     void setSites(const arma::mat &points);
 
     // Set the value to be colorScaleped in each site
     void setValues(const arma::vec &values);
 
-    // Set colorScale data based on the given color scale
-    void setColorScale(std::shared_ptr<const ColorScale> scale);
+    // Set colormap texture used for color lookup in rendering
+    void setColormap(GLuint texture);
 
     void setScale(const LinearScale<float> &sx,
                   const LinearScale<float> &sy);
@@ -59,11 +71,30 @@ public:
     // Maximum blur radius
     void setBeta(float beta);
 
+    void setupShaders();
+    void setupVAOs();
+    void setupTextures();
+    void resizeTextures();
+
+    void updateSites();
+    void updateValues();
+    void updateColormap();
+    void updateTransform();
+
 private:
-    std::vector<float> m_sites, m_values, m_cmap;
-    LinearScale<float> m_sx, m_sy;
+    std::unique_ptr<Shader> m_program1, m_program2, m_programVoronoi;
+    GLuint m_dtVAO, m_sitesVAO, m_2ndPassVAO, m_voronoiVAO;
+    GLfloat m_transform[4][4];
+    GLuint m_FBO, m_voronoiFBO, m_preFBO;
+    GLuint m_VBOs[3], m_dtVBO;
+    GLuint m_textures[2], m_colormapTex, m_voronoiDepthTex, m_outTex;
+
     float m_alpha, m_beta;
-    bool m_sitesChanged, m_valuesChanged, m_colorScaleChanged;
+    size_t m_width, m_height;
+    std::vector<vec2> m_sites;
+    std::vector<float> m_values, m_cmap;
+    LinearScale<float> m_sx, m_sy;
+    bool m_sitesChanged, m_valuesChanged, m_colorScaleChanged, m_redraw;
 };
 
 #endif // VORONOISPLAT_H