diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3f0f921..50375e4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -48,3 +48,13 @@ target_link_libraries(fmri ${png++_LIBRARIES})
 find_package(GLUT REQUIRED)
 include_directories(${GLUT_INCLUDE_DIRS})
 target_link_libraries(fmri ${GLUT_LIBRARIES})
+
+# Require OpenGL (visualisation)
+find_package(OpenGL REQUIRED)
+include_directories(${OPENGL_INCLUDE_DIRS})
+target_link_libraries(fmri ${OPENGL_LIBRARIES})
+
+# Require GLEW (visulalisation)
+find_package(GLEW REQUIRED)
+include_directories(${GLEW_INCLUDE_DIRS})
+target_link_libraries(fmri ${GLEW_LIBRARIES})
diff --git a/src/glutils.cpp b/src/glutils.cpp
new file mode 100644
index 0000000..02d0b65
--- /dev/null
+++ b/src/glutils.cpp
@@ -0,0 +1,104 @@
+#include <cstdint>
+#include <memory>
+#include <vector>
+#include <GL/glew.h>
+#include <cstring>
+#include <glog/logging.h>
+#include "glutils.hpp"
+
+using namespace fmri;
+using namespace std;
+
+/**
+ * Utility function to send OpenGL info dumps to the log
+ *
+ * @param object The GLObject to get info for
+ * @param glGet__iv The function to query log length with
+ * @param glGet__InfoLog The function to get the log with
+ */
+static void show_info_log(
+        GLuint object,
+        PFNGLGETSHADERIVPROC glGet__iv,
+        PFNGLGETSHADERINFOLOGPROC glGet__InfoLog
+)
+{
+    GLint log_length;
+    glGet__iv(object, GL_INFO_LOG_LENGTH, &log_length);
+    unique_ptr<char[]> log(new char[log_length]);
+    glGet__InfoLog(object, log_length, nullptr, log.get());
+
+    LOG(INFO) << log.get() << endl;
+}
+
+GLuint fmri::loadTexture(DType const *data, unsigned int width, unsigned int height)
+{
+    // Load and scale texture
+    vector<DType> textureBuffer(data, data + (width * height));
+    rescale(textureBuffer.begin(), textureBuffer.end(), 0, 1);
+
+    const float color[] = {0, 0, 0}; // Background color for textures.
+    GLuint texture;
+    glGenTextures(1, &texture);
+
+    glBindTexture(GL_TEXTURE_2D, texture);
+
+    // Set up (lack of) repetition
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
+    glTexParameterfv(GL_TEXTURE_2D, GL_TEXTURE_BORDER_COLOR, color);
+
+    static_assert(sizeof(DType) == 4); // verify data type for texture.
+    glTexImage2D(GL_TEXTURE_2D, 0, GL_R32F, 2, 2, 0, GL_R32F, GL_FLOAT, textureBuffer.data());
+
+    // Set up texture scaling
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR); // Use mipmapping for scaling down
+    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST); // Use nearest pixel when scaling up.
+
+    glGenerateMipmap(GL_TEXTURE_2D);
+
+
+    return texture;
+}
+
+void fmri::changeWindowSize(int w, int h)
+{
+    // Prevent a divide by zero, when window is too short
+    // (you cant make a window of zero width).
+    if (h == 0)
+        h = 1;
+
+    float ratio = w * 1.0f / h;
+
+    // Use the Projection Matrix
+    glMatrixMode(GL_PROJECTION);
+
+    // Reset Matrix
+    glLoadIdentity();
+
+    // Set the viewport to be the entire window
+    glViewport(0, 0, w, h);
+
+    // Set the correct perspective.
+    gluPerspective(45.0f, ratio, 0.1f, 100.0f);
+
+    // Get Back to the Modelview
+    glMatrixMode(GL_MODELVIEW);
+}
+
+GLuint fmri::compileShader(GLenum type, char const *source)
+{
+    GLuint shader = glCreateShader(type);
+    auto length = static_cast<GLint>(strlen(source));
+    GLint compileOK;
+
+    glShaderSource(shader, 1, &source, &length);
+    glCompileShader(shader);
+
+    glGetShaderiv(shader, GL_COMPILE_STATUS, &compileOK);
+    if (!compileOK) {
+        LOG(WARNING) << "Failed to compile " << source << endl;
+        show_info_log(shader, glGetShaderiv, glGetShaderInfoLog);
+        glDeleteShader(shader);
+        abort();
+    }
+}
diff --git a/src/glutils.hpp b/src/glutils.hpp
new file mode 100644
index 0000000..6ee8b05
--- /dev/null
+++ b/src/glutils.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+#include "LayerData.hpp"
+#include "utils.hpp"
+#include <GL/glut.h>
+
+namespace fmri {
+    /**
+     * Create a texture from an array of data.
+     *
+     * @param data
+     * @param width
+     * @param height
+     * @return A texture reference.
+     */
+    GLuint loadTexture(DType const * data, unsigned int width, unsigned int height);
+
+    /**
+     *
+     * @param type
+     * @param source
+     * @return
+     */
+    GLuint compileShader(GLenum type, char const * source);
+
+    /**
+     * Callback handler to handle resizing windows.
+     *
+     * This function resizes the rendering viewport so everything still
+     * looks proportional.
+     *
+     * @param w new Width
+     * @param h new Height.
+     */
+    void changeWindowSize(int w, int h);
+}
diff --git a/src/main.cpp b/src/main.cpp
index a147171..8a7aada 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,44 +1,135 @@
 #include <algorithm>
 #include <iostream>
 #include <glog/logging.h>
+#include <GL/glew.h>
+#include <GL/glut.h>
+#include <map>
+#include <tuple>
 
+#include "LayerData.hpp"
 #include "Options.hpp"
 #include "Simulator.hpp"
+#include "glutils.hpp"
 
 using namespace std;
 using namespace fmri;
 
-int main(int argc, char *const argv[]) {
-    google::InitGoogleLogging(argv[0]);
+struct {
+    optional<vector<string>> labels;
+    vector<vector<LayerData>> data;
+    float angle = 0;
+    map<tuple<string, int, int>, GLuint> textureMap;
+} rendererData;
 
-    Options options = Options::parse(argc, argv);
-    auto labels = options.labels();
+static vector<vector<LayerData>> getSimulationData(const Options& options)
+{
+    vector<vector<LayerData>> results;
     auto dumper = options.imageDumper();
-
     Simulator simulator(options.model(), options.weights(), options.means());
 
-	for (const auto &image : options.inputs()) {
-		const auto res = simulator.simulate(image);
-		LOG(INFO) << "Result for " << image << ":" << endl;
+    for (const auto& image : options.inputs()) {
+        results.emplace_back(simulator.simulate(image));
+    }
 
-		const auto& resultRow = res[res.size() - 1];
-		if (labels) {
-			vector<DType> weights(resultRow.data(), resultRow.data() + resultRow.numEntries());
-			auto scores = combine(weights, *labels);
-			sort(scores.begin(), scores.end(), greater<>());
-			for (unsigned int i = 0; i < scores.size() && i < 5; ++i) {
-				LOG(INFO) << scores[i].first << " " << scores[i].second << endl;
-			}
-		} else {
-			LOG(INFO) << "Best result: " << *(resultRow.data(), resultRow.data() + resultRow.numEntries()) << endl;
-		}
+    CHECK_GT(results.size(), 0) << "Should have some results" << endl;
 
-        if (dumper) {
-            for (const auto& layer : res) {
-                dumper->dump(layer);
+    if (dumper) {
+        for (auto& layer : *results.begin()) {
+            dumper->dump(layer);
+        }
+    }
+
+    return results;
+}
+
+static void render()
+{
+
+    // Clear Color and Depth Buffers
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+    // Reset transformations
+    glLoadIdentity();
+    // Set the camera
+    gluLookAt(	0.0f, 0.0f, 10.0f,
+                  0.0f, 0.0f,  0.0f,
+                  0.0f, 1.0f,  0.0f);
+
+    glRotatef(rendererData.angle, 0.0f, 1.0f, 0.0f);
+
+    glBegin(GL_TRIANGLES);
+    glVertex3f(-2.0f,-2.0f, 0.0f);
+    glVertex3f( 2.0f, 0.0f, 0.0);
+    glVertex3f( 0.0f, 2.0f, 0.0);
+    glEnd();
+
+    rendererData.angle+=0.1f;
+
+    glutSwapBuffers();
+}
+
+static void reloadTextures(unsigned dataIndex)
+{
+    // First, release any existing textures
+    for (auto& entry : rendererData.textureMap) {
+        glDeleteTextures(0, &entry.second);
+    }
+
+    rendererData.textureMap.clear();
+
+    const auto& dataSet = rendererData.data[dataIndex];
+
+    for (auto& layer : dataSet) {
+        auto dimensions = layer.shape();
+        if (dimensions.size() != 4) {
+            continue;
+        }
+
+        const auto images = dimensions[0],
+                channels = dimensions[1],
+                width = dimensions[2],
+                height = dimensions[3];
+
+        auto dataPtr = layer.data();
+        for (auto i = 0; i < images; ++i) {
+            for (auto j = 0; j < channels; ++j) {
+                rendererData.textureMap[make_tuple(layer.name(), i, j)] = loadTexture(dataPtr, width, height);
+                dataPtr += width * height;
             }
         }
-	}
+    }
+}
+
+
+
+int main(int argc, char * argv[]) {
+    google::InitGoogleLogging(argv[0]);
+    google::InstallFailureSignalHandler();
+
+    glutInit(&argc, argv);
+    glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA);
+    glutCreateWindow(argv[0]);
+
+    // Prepare data for simulations
+    Options options = Options::parse(argc, argv);
+    rendererData.labels =  options.labels();
+    rendererData.data = getSimulationData(options);
+
+    // Register callbacks
+    glutDisplayFunc(render);
+    glutIdleFunc(render);
+    glutReshapeFunc(changeWindowSize);
+
+    glewInit();
+    if (!GLEW_VERSION_2_0) {
+        cerr << "OpenGL 2.0 not available" << endl;
+        return 2;
+    }
+
+    reloadTextures(0);
+
+    // Start visualisation
+    glutMainLoop();
 
     google::ShutdownGoogleLogging();
 
diff --git a/src/shaders.cpp b/src/shaders.cpp
new file mode 100644
index 0000000..2ded97d
--- /dev/null
+++ b/src/shaders.cpp
@@ -0,0 +1 @@
+#include "shaders.hpp"
diff --git a/src/shaders.hpp b/src/shaders.hpp
new file mode 100644
index 0000000..a34892a
--- /dev/null
+++ b/src/shaders.hpp
@@ -0,0 +1,5 @@
+#pragma once
+
+#include <GL/glut.h>
+
+extern const char * GRAYSCALE_TEXTURE_SHADER;