Restructure project sources.

src/fmri/ActivityAnimation.cpp

@@ -0,0 +1,68 @@
+#include <cmath>
+#include <cstring>
+#include <random>
+#include <GL/gl.h>
+#include <caffe/util/math_functions.hpp>
+#include "Range.hpp"
+#include "ActivityAnimation.hpp"
+
+using namespace std;
+using namespace fmri;
+
+ActivityAnimation::Color ActivityAnimation::colorBySign(float intensity)
+{
+    if (intensity > 0) {
+        return {0, 1, 0};
+    } else {
+        return {1, 0, 0};
+    }
+}
+
+ActivityAnimation::ActivityAnimation(
+            const std::vector<std::pair<DType, std::pair<std::size_t, std::size_t>>> &interactions,
+            const float *aPositions, const float *bPositions) :
+        ActivityAnimation(interactions, aPositions, bPositions, ActivityAnimation::colorBySign)
+{
+}
+
+ActivityAnimation::ActivityAnimation(
+            const std::vector<std::pair<DType, std::pair<std::size_t, std::size_t>>> &interactions,
+            const float *aPositions, const float *bPositions, ColoringFunction coloring)
+        :
+        bufferLength(3 * interactions.size())
+{
+    CHECK(coloring) << "Invalid coloring function passed.";
+    startingPos.reserve(bufferLength);
+    delta.reserve(bufferLength);
+    colorBuf.reserve(interactions.size());
+
+    for (auto &entry : interactions) {
+        auto *aPos = &aPositions[3 * entry.second.first];
+        auto *bPos = &bPositions[3 * entry.second.second];
+
+        colorBuf.push_back(coloring(entry.first));
+
+        for (auto i : Range(3)) {
+            startingPos.emplace_back(aPos[i]);
+            delta.emplace_back(bPos[i] - aPos[i] + (i % 3 ? 0 : LAYER_X_OFFSET));
+        }
+    }
+}
+
+void ActivityAnimation::draw(float timeScale)
+{
+    std::unique_ptr<float[]> vertexBuffer(new float[bufferLength]);
+    caffe::caffe_copy(bufferLength, delta.data(), vertexBuffer.get());
+    caffe::caffe_scal(bufferLength, timeScale, vertexBuffer.get());
+    caffe::caffe_add(bufferLength, startingPos.data(), vertexBuffer.get(), vertexBuffer.get());
+
+    glPointSize(5);
+
+    glEnableClientState(GL_VERTEX_ARRAY);
+    glEnableClientState(GL_COLOR_ARRAY);
+    glColorPointer(3, GL_FLOAT, 0, colorBuf.data());
+    glVertexPointer(3, GL_FLOAT, 0, vertexBuffer.get());
+    glDrawArrays(GL_POINTS, 0, bufferLength / 3);
+    glDisableClientState(GL_COLOR_ARRAY);
+    glDisableClientState(GL_VERTEX_ARRAY);
+}

src/fmri/ActivityAnimation.hpp

@@ -0,0 +1,34 @@
+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <vector>
+#include "Animation.hpp"
+#include "utils.hpp"
+
+namespace fmri
+{
+    class ActivityAnimation
+            : public Animation
+    {
+    public:
+        typedef std::array<float, 3> Color;
+        typedef std::function<Color(float)> ColoringFunction;
+
+        ActivityAnimation(
+                    const std::vector<std::pair<DType, std::pair<std::size_t, std::size_t>>> &interactions,
+                    const float *aPositions, const float *bPositions);
+        ActivityAnimation(
+                    const std::vector<std::pair<DType, std::pair<std::size_t, std::size_t>>> &interactions,
+                    const float *aPositions, const float *bPositions, ColoringFunction coloring);
+        void draw(float timeScale) override;
+
+        static Color colorBySign(float intensity);
+
+    private:
+        std::size_t bufferLength;
+        std::vector<std::array<float, 3>> colorBuf;
+        std::vector<float> startingPos;
+        std::vector<float> delta;
+    };
+}

src/fmri/Animation.cpp

@@ -0,0 +1,5 @@
+//
+// Created by bert on 13/02/18.
+//
+
+#include "Animation.hpp"

src/fmri/Animation.hpp

@@ -0,0 +1,14 @@
+#pragma once
+
+
+namespace fmri
+{
+    class Animation
+    {
+    public:
+        virtual ~Animation() = default;
+
+        virtual void draw(float step) = 0;
+
+    };
+}

src/fmri/DummyLayerVisualisation.hpp

@@ -0,0 +1,16 @@
+#pragma once
+
+#include "LayerVisualisation.hpp"
+
+namespace fmri
+{
+    /**
+     * Visualisation that does not actually do anything.
+     */
+    class DummyLayerVisualisation : public LayerVisualisation
+    {
+    public:
+        void render() override
+        {};
+    };
+}

src/fmri/FlatLayerVisualisation.cpp

@@ -0,0 +1,101 @@
+#include <glog/logging.h>
+#include <GL/gl.h>
+
+#include "FlatLayerVisualisation.hpp"
+#include "Range.hpp"
+
+using namespace std;
+using namespace fmri;
+
+static inline void computeColor(float intensity, float limit, float *destination)
+{
+    const float saturation = min(-log(abs(intensity) / limit) / 10.0f, 1.0f);
+    if (intensity > 0) {
+        destination[0] = saturation;
+        destination[1] = saturation;
+        destination[2] = 1;
+    } else {
+        destination[0] = 1;
+        destination[1] = saturation;
+        destination[2] = saturation;
+    }
+}
+
+FlatLayerVisualisation::FlatLayerVisualisation(const LayerData &layer, Ordering ordering) :
+        LayerVisualisation(layer.numEntries()),
+        ordering(ordering),
+        faceCount(layer.numEntries() * NODE_FACES.size() / 3),
+        vertexBuffer(new float[faceCount * 3]),
+        colorBuffer(new float[faceCount * 3]),
+        indexBuffer(new int[faceCount * 3])
+{
+    auto &shape = layer.shape();
+    CHECK_EQ(shape.size(), 2) << "layer should be flat!" << endl;
+    CHECK_EQ(shape[0], 1) << "Only single images supported." << endl;
+
+    initializeNodePositions();
+
+    const auto limit = (int) layer.numEntries();
+    auto data = layer.data();
+    const auto
+    [minElem, maxElem] = minmax_element(data, data + limit);
+
+    auto scalingMax = max(abs(*minElem), abs(*maxElem));
+
+    int v = 0;
+    for (int i : Range(limit)) {
+        setVertexPositions(i, vertexBuffer.get() + NODE_FACES.size() * i);
+        const auto vertexBase = static_cast<int>(i * NODE_FACES.size() / 3);
+
+        // Define the colors for the vertices
+        for (auto c : Range(NODE_SHAPE.size() / 3)) {
+            computeColor(data[i], scalingMax, &colorBuffer[NODE_FACES.size() * i + 3 * c]);
+        }
+
+        // Set the face nodes indices
+        for (auto faceNode : NODE_FACES) {
+            indexBuffer[v++] = vertexBase + faceNode;
+        }
+    }
+    assert(v == (int) faceCount * 3);
+}
+
+void FlatLayerVisualisation::render()
+{
+    glEnableClientState(GL_VERTEX_ARRAY);
+    glEnableClientState(GL_COLOR_ARRAY);
+
+    glVertexPointer(3, GL_FLOAT, 0, vertexBuffer.get());
+    glColorPointer(3, GL_FLOAT, 0, colorBuffer.get());
+    glDrawElements(GL_TRIANGLES, faceCount * 3, GL_UNSIGNED_INT, indexBuffer.get());
+    glDisableClientState(GL_COLOR_ARRAY);
+
+    // Now draw wireframe
+    glColor4f(0, 0, 0, 1);
+    glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
+    glDrawElements(GL_TRIANGLES, faceCount * 3, GL_UNSIGNED_INT, indexBuffer.get());
+    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
+
+
+    glDisableClientState(GL_VERTEX_ARRAY);
+}
+
+void FlatLayerVisualisation::setVertexPositions(const int vertexNo, float *destination)
+{
+    for (auto i : Range(NODE_SHAPE.size())) {
+        destination[i] = NODE_SHAPE[i] + nodePositions_[3 * vertexNo + (i % 3)];
+    }
+}
+
+void FlatLayerVisualisation::initializeNodePositions()
+{
+    switch (ordering) {
+        case Ordering::LINE:
+            initNodePositions<Ordering::LINE>(faceCount / 4, 2);
+            break;
+
+        case Ordering::SQUARE:
+            initNodePositions<Ordering::SQUARE>(faceCount / 4, 2);
+            break;
+    }
+}

src/fmri/FlatLayerVisualisation.hpp

@@ -0,0 +1,42 @@
+#pragma once
+
+#include <memory>
+
+#include "LayerData.hpp"
+#include "LayerVisualisation.hpp"
+
+namespace fmri
+{
+    class FlatLayerVisualisation : public LayerVisualisation
+    {
+    public:
+        explicit FlatLayerVisualisation(const LayerData &layer, Ordering ordering);
+
+        void render() override;
+
+    private:
+        Ordering ordering;
+        std::size_t faceCount;
+        std::unique_ptr<float[]> vertexBuffer;
+        std::unique_ptr<float[]> colorBuffer;
+        std::unique_ptr<int[]> indexBuffer;
+
+        static constexpr const std::array<float, 12> NODE_SHAPE = {
+                -0.5f, 0, 0.5f,
+                0, 0, -0.5f,
+                0, 1, 0,
+                0.5f, 0, 0.5f
+        };
+        static constexpr const std::array<int, 12> NODE_FACES = {
+                0, 1, 2,
+                0, 1, 3,
+                0, 2, 3,
+                1, 2, 3
+        };
+
+        void setVertexPositions(int vertexNo, float *destination);
+
+        // Various functions defining the way the nodes will be aligned.
+        void initializeNodePositions();
+    };
+}

src/fmri/ImageInteractionAnimation.cpp

@@ -0,0 +1,30 @@
+#include "ImageInteractionAnimation.hpp"
+#include "glutils.hpp"
+#include "MultiImageVisualisation.hpp"
+#include <caffe/util/math_functions.hpp>
+
+using namespace fmri;
+
+
+void ImageInteractionAnimation::draw(float step)
+{
+    auto vertexBuffer = deltas;
+    caffe::caffe_scal(deltas.size(), step, vertexBuffer.data());
+    caffe::caffe_add(vertexBuffer.size(), vertexBuffer.data(), startingPositions.data(), vertexBuffer.data());
+
+    drawImageTiles(vertexBuffer.size() / 3, vertexBuffer.data(), textureCoordinates.data(), texture);
+}
+
+ImageInteractionAnimation::ImageInteractionAnimation(const DType *data, const std::vector<int> &shape, const std::vector<float> &prevPositions,
+                                                     const std::vector<float> &curPositions) :
+        texture(loadTexture(data, shape[2], shape[1] * shape[3], shape[1])),
+        startingPositions(MultiImageVisualisation::getVertices(prevPositions)),
+        deltas(MultiImageVisualisation::getVertices(curPositions)),
+        textureCoordinates(MultiImageVisualisation::getTexCoords(shape[1]))
+{
+    caffe::caffe_sub(deltas.size(), deltas.data(), startingPositions.data(), deltas.data());
+
+    for (auto i = 0u; i < deltas.size(); i += 3) {
+        deltas[i] = LAYER_X_OFFSET;
+    }
+}

src/fmri/ImageInteractionAnimation.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include "Animation.hpp"
+#include "utils.hpp"
+#include "Texture.hpp"
+
+namespace fmri
+{
+    class ImageInteractionAnimation : public Animation
+    {
+    public:
+        ImageInteractionAnimation(const DType *data, const std::vector<int> &shape, const std::vector<float> &prevPositions,
+                                          const std::vector<float> &curPositions);
+        virtual void draw(float step);
+
+    private:
+        Texture texture;
+        std::vector<float> startingPositions;
+        std::vector<float> deltas;
+        std::vector<float> textureCoordinates;
+    };
+}

src/fmri/InputLayerVisualisation.cpp

@@ -0,0 +1,95 @@
+#include <caffe/util/math_functions.hpp>
+#include <GL/glu.h>
+#include <opencv2/core/mat.hpp>
+#include <opencv2/core.hpp>
+#include "InputLayerVisualisation.hpp"
+#include "Range.hpp"
+
+using namespace fmri;
+using namespace std;
+
+/**
+ * Combine an arbitrary number of channels into an RGB image.
+ *
+ * If there are less than 3 channels, the first channel is repeated. This
+ * results in greyscale images for single-channel images. Any channels after
+ * the third are ignored.
+ *
+ * @param data Layer data to generate an image for.
+ * @return A normalized RGB image, stored in a vector.
+ */
+static vector<float> getRGBImage(const LayerData &data)
+{
+    vector<cv::Mat> channels;
+    const int numPixels = data.shape()[2] * data.shape()[3];
+    for (auto i : Range(3)) {
+        if (i >= data.shape()[1]) {
+            channels.push_back(channels[0]);
+        }
+
+        cv::Mat channel(data.shape()[3], data.shape()[2], CV_32FC1);
+        copy(data.data() + i * numPixels, data.data() + (i + 1) * numPixels, channel.begin<float>());
+        channels.push_back(channel);
+    }
+
+    swap(channels[0], channels[2]);
+
+    cv::Mat outImage;
+    cv::merge(channels, outImage);
+
+    outImage = outImage.reshape(1);
+
+    vector<float> final(outImage.begin<float>(), outImage.end<float>());
+    rescale(final.begin(), final.end(), 0.f, 1.f);
+
+    return final;
+}
+
+InputLayerVisualisation::InputLayerVisualisation(const LayerData &data)
+{
+    CHECK_EQ(data.shape().size(), 4) << "Should be image-like-layer." << endl;
+    auto imageData = getRGBImage(data);
+
+    const auto images = data.shape()[0], channels = data.shape()[1], width = data.shape()[2], height = data.shape()[3];
+    CHECK_EQ(images, 1) << "Should be single image" << endl;
+
+    targetWidth = width / 5.f;
+    targetHeight = width / 5.f;
+
+    nodePositions_ = {0, targetHeight / 2, targetWidth / -2};
+    for (auto i : Range(3, 3 * channels)) {
+        nodePositions_.push_back(nodePositions_[i % 3]);
+    }
+
+    texture.configure(GL_TEXTURE_2D);
+    gluBuild2DMipmaps(GL_TEXTURE_2D, GL_RGB, width, height, GL_RGB, GL_FLOAT, imageData.data());
+}
+
+void InputLayerVisualisation::render()
+{
+    const float vertices[] = {
+            0, 0, 0,
+            0, 0, -targetWidth,
+            0, targetHeight, -targetWidth,
+            0, targetHeight, 0,
+    };
+
+    const float texCoords[] = {
+            0, 1,
+            1, 1,
+            1, 0,
+            0, 0,
+    };
+
+    glEnableClientState(GL_VERTEX_ARRAY);
+    glEnableClientState(GL_TEXTURE_COORD_ARRAY);
+    glVertexPointer(3, GL_FLOAT, 0, vertices);
+    glTexCoordPointer(2, GL_FLOAT, 0, texCoords);
+    glEnable(GL_TEXTURE_2D);
+    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
+    texture.bind(GL_TEXTURE_2D);
+    glDrawArrays(GL_QUADS, 0, 4);
+    glDisable(GL_TEXTURE_2D);
+    glDisableClientState(GL_TEXTURE_COORD_ARRAY);
+    glDisableClientState(GL_VERTEX_ARRAY);
+}

src/fmri/InputLayerVisualisation.hpp

@@ -0,0 +1,22 @@
+#pragma once
+
+#include "LayerData.hpp"
+#include "LayerVisualisation.hpp"
+#include "Texture.hpp"
+
+namespace fmri
+{
+    class InputLayerVisualisation : public LayerVisualisation
+    {
+    public:
+        explicit InputLayerVisualisation(const LayerData &data);
+
+        void render() override;
+
+    private:
+        Texture texture;
+        float targetWidth;
+        float targetHeight;
+
+    };
+}

src/fmri/LayerData.cpp

@@ -0,0 +1,63 @@
+#include <cstring>
+#include <functional>
+#include <iostream>
+#include <numeric>
+
+#include <glog/logging.h>
+
+#include "LayerData.hpp"
+
+using namespace fmri;
+using namespace std;
+
+LayerData::LayerData(const string& name, const vector<int>& shape, const DType* data) :
+	name_(name),
+	shape_(shape)
+{
+	const auto dataSize = numEntries();
+	// Compute the dimension of the data area
+	data_.reset(new DType[dataSize]);
+
+	// Copy the data over with memcpy because it's just faster that way
+	memcpy(data_.get(), data, sizeof(DType) * dataSize);
+}
+
+size_t LayerData::numEntries() const
+{
+	return static_cast<size_t>(accumulate(shape_.begin(), shape_.end(), 1, multiplies<>()));
+}
+
+const vector<int>& LayerData::shape() const
+{
+	return shape_;
+}
+
+const string& LayerData::name() const
+{
+	return name_;
+}
+
+DType const * LayerData::data() const
+{
+	return data_.get();
+}
+
+ostream& operator<< (ostream& o, const LayerData& layer)
+{
+    o << layer.name() << '(';
+    bool first = true;
+
+    for (auto d : layer.shape()) {
+        if (!first) {
+            o << ", ";
+        } else {
+            first = false;
+        }
+
+        o << d;
+    }
+
+    o << ')';
+
+    return o;
+}

src/fmri/LayerData.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <iostream>
+#include <memory>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "utils.hpp"
+
+namespace fmri
+{
+
+    using std::ostream;
+    using std::string;
+    using std::string_view;
+    using std::unique_ptr;
+    using std::vector;
+
+    class LayerData
+    {
+    public:
+        LayerData(const string &name, const vector<int> &shape, const DType *data);
+        LayerData(const LayerData &) = delete;
+
+        LayerData(LayerData &&) = default;
+        LayerData &operator=(const LayerData &) = delete;
+        LayerData &operator=(LayerData &&) = default;
+
+        const string &name() const;
+        const vector<int> &shape() const;
+        DType const *data() const;
+        size_t numEntries() const;
+    private:
+        string name_;
+        vector<int> shape_;
+        unique_ptr<DType[]> data_;
+    };
+}
+
+std::ostream& operator<<(std::ostream&, const fmri::LayerData&);

src/fmri/LayerInfo.cpp

@@ -0,0 +1,62 @@
+#include "LayerInfo.hpp"
+
+using namespace std;
+using namespace fmri;
+
+const unordered_map<string_view, LayerInfo::Type> LayerInfo::NAME_TYPE_MAP = {
+        {"Input", Type::Input},
+        {"Convolution", Type::Convolutional},
+        {"ReLU", Type::ReLU},
+        {"Pooling", Type::Pooling},
+        {"InnerProduct", Type::InnerProduct},
+        {"Dropout", Type::DropOut},
+        {"LRN", Type::LRN},
+        {"Split", Type::Split},
+        {"Softmax", Type::Softmax}
+};
+
+
+LayerInfo::Type LayerInfo::typeByName(string_view name)
+{
+    try {
+        return NAME_TYPE_MAP.at(name);
+    } catch (std::out_of_range &e) {
+        LOG(INFO) << "Received unknown layer type: " << name << endl;
+        return Type::Other;
+    }
+}
+
+LayerInfo::LayerInfo(string_view name, string_view type,
+                     const vector<boost::shared_ptr<caffe::Blob<DType>>> &parameters)
+: parameters_(parameters), type_(typeByName(type)), name_(name)
+{
+}
+
+const std::string &LayerInfo::name() const
+{
+    return name_;
+}
+
+LayerInfo::Type LayerInfo::type() const
+{
+    return type_;
+}
+
+const std::vector<boost::shared_ptr<caffe::Blob<DType>>>& LayerInfo::parameters() const
+{
+    return parameters_;
+}
+
+std::ostream &fmri::operator<<(std::ostream &out, LayerInfo::Type type)
+{
+    for (auto i : LayerInfo::NAME_TYPE_MAP) {
+        if (i.second == type) {
+            out << i.first;
+            return out;
+        }
+    }
+
+    out << "ERROR! UNSUPPORTED TYPE";
+    return out;
+}
+

src/fmri/LayerInfo.hpp

@@ -0,0 +1,47 @@
+#pragma once
+
+#include <string_view>
+#include <caffe/blob.hpp>
+#include <string>
+#include "utils.hpp"
+
+namespace fmri
+{
+    class LayerInfo
+    {
+    public:
+        enum class Type
+        {
+            Input,
+            Convolutional,
+            ReLU,
+            Pooling,
+            InnerProduct,
+            DropOut,
+            LRN,
+            Split,
+            Softmax,
+            Other
+        };
+
+        LayerInfo(std::string_view name, std::string_view type,
+                  const std::vector<boost::shared_ptr<caffe::Blob<DType>>> &parameters);
+
+        const std::string& name() const;
+        Type type() const;
+        const std::vector<boost::shared_ptr<caffe::Blob<DType>>>& parameters() const;
+
+        static Type typeByName(std::string_view name);
+
+        friend std::ostream& operator<<(std::ostream& out, Type type);
+
+    private:
+        std::vector<boost::shared_ptr<caffe::Blob<DType>>> parameters_;
+        Type type_;
+        std::string name_;
+
+        const static std::unordered_map<std::string_view, Type> NAME_TYPE_MAP;
+    };
+
+    std::ostream& operator<<(std::ostream& out, LayerInfo::Type type);
+}

src/fmri/LayerVisualisation.cpp

@@ -0,0 +1,40 @@
+#include "LayerVisualisation.hpp"
+#include "Range.hpp"
+#include "utils.hpp"
+
+const std::vector<float> &fmri::LayerVisualisation::nodePositions() const
+{
+    return nodePositions_;
+}
+
+fmri::LayerVisualisation::LayerVisualisation(size_t numNodes)
+        : nodePositions_(numNodes * 3)
+{
+}
+
+template<>
+void fmri::LayerVisualisation::initNodePositions<fmri::LayerVisualisation::Ordering::LINE>(size_t n, float spacing)
+{
+    nodePositions_.clear();
+    nodePositions_.reserve(3 * n);
+
+    for (auto i : Range(n)) {
+        nodePositions_.push_back(0);
+        nodePositions_.push_back(0);
+        nodePositions_.push_back(-spacing * i);
+    }
+}
+
+template<>
+void fmri::LayerVisualisation::initNodePositions<fmri::LayerVisualisation::Ordering::SQUARE>(size_t n, float spacing)
+{
+    nodePositions_.clear();
+    nodePositions_.reserve(3 * n);
+    const auto columns = numCols(n);
+
+    for (auto i : Range(n)) {
+        nodePositions_.push_back(0);
+        nodePositions_.push_back(spacing * (i / columns));
+        nodePositions_.push_back(-spacing * (i % columns));
+    }
+}

src/fmri/LayerVisualisation.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <vector>
+
+namespace fmri
+{
+    class LayerVisualisation
+    {
+    public:
+        enum class Ordering {
+            LINE,
+            SQUARE,
+        };
+
+        LayerVisualisation() = default;
+        explicit LayerVisualisation(size_t numNodes);
+        virtual ~LayerVisualisation() = default;
+
+        virtual void render() = 0;
+        virtual const std::vector<float>& nodePositions() const;
+
+    protected:
+        std::vector<float> nodePositions_;
+
+        template<Ordering Order>
+        void initNodePositions(size_t n, float spacing);
+    };
+}

src/fmri/MultiImageVisualisation.cpp

@@ -0,0 +1,67 @@
+#include <glog/logging.h>
+#include "MultiImageVisualisation.hpp"
+#include "glutils.hpp"
+#include "Range.hpp"
+
+using namespace fmri;
+using namespace std;
+
+MultiImageVisualisation::MultiImageVisualisation(const fmri::LayerData &layer)
+{
+    auto dimensions = layer.shape();
+    CHECK_EQ(4, dimensions.size()) << "Should be image-like layer";
+
+    const auto images = dimensions[0],
+            channels = dimensions[1],
+            width = dimensions[2],
+            height = dimensions[3];
+
+    CHECK_EQ(1, images) << "Only single input image is supported" << endl;
+
+    texture = loadTexture(layer.data(), width, channels * height, channels);
+    initNodePositions<Ordering::SQUARE>(channels, 3);
+    vertexBuffer = getVertices(nodePositions_);
+    texCoordBuffer = getTexCoords(channels);
+}
+
+void MultiImageVisualisation::render()
+{
+    drawImageTiles(vertexBuffer.size() / 3, vertexBuffer.data(), texCoordBuffer.data(), texture);
+}
+
+vector<float> MultiImageVisualisation::getVertices(const std::vector<float> &nodePositions, float scaling)
+{
+    std::vector<float> vertices;
+    vertices.reserve(nodePositions.size() * BASE_VERTICES.size() / 3);
+    for (auto i = 0u; i < nodePositions.size(); i += 3) {
+        auto pos = &nodePositions[i];
+        for (auto j = 0u; j < BASE_VERTICES.size(); ++j) {
+            vertices.push_back(BASE_VERTICES[j] * scaling + pos[j % 3]);
+        }
+    }
+
+    return vertices;
+}
+
+std::vector<float> MultiImageVisualisation::getTexCoords(int n)
+{
+    std::vector<float> coords;
+    coords.reserve(8 * n);
+
+    const float channels = n;
+
+    for (int i = 0; i < n; ++i) {
+        std::array<float, 8> textureCoords = {
+                1, (i + 1) / channels,
+                1, i / channels,
+                0, i / channels,
+                0, (i + 1) / channels,
+        };
+
+        for (auto coord : textureCoords) {
+            coords.push_back(coord);
+        }
+    }
+
+    return coords;
+}

src/fmri/MultiImageVisualisation.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include <map>
+#include <memory>
+#include "LayerVisualisation.hpp"
+#include "LayerData.hpp"
+#include "Texture.hpp"
+
+namespace fmri
+{
+    class MultiImageVisualisation : public LayerVisualisation
+    {
+    public:
+        constexpr const static std::array<float, 12> BASE_VERTICES = {
+                0, -1, -1,
+                0, 1, -1,
+                0, 1, 1,
+                0, -1, 1,
+        };
+
+        explicit MultiImageVisualisation(const LayerData&);
+
+        void render() override;
+
+        static vector<float> getVertices(const std::vector<float> &nodePositions, float scaling = 1);
+        static std::vector<float> getTexCoords(int n);
+
+    private:
+        Texture texture;
+        std::vector<float> vertexBuffer;
+        std::vector<float> texCoordBuffer;
+    };
+}

src/fmri/Options.cpp

@@ -0,0 +1,147 @@
+#include <algorithm>
+#include <cstdio>
+#include <iostream>
+#include <unistd.h>
+#include "Options.hpp"
+#include "utils.hpp"
+
+using namespace fmri;
+using namespace std;
+
+static void show_help(const char *progname, int exitcode) {
+    cerr << "Usage: " << progname << " -m MODEL -w WEIGHTS INPUTS..." << endl
+         << endl
+         << R"END(Simulate the specified network on the specified inputs.
+
+Options:
+	-h	show this message
+	-n	(required) the model file to simulate
+	-w	(required) the trained weights
+    -m  means file. Will be substracted from input if available.
+    -l  labels file. Will be used to print prediction labels if available.
+    -d  Image dump dir. Will be filled with PNG images of intermediate results.)END" << endl;
+
+    exit(exitcode);
+}
+
+static void check_file(const char *filename) {
+    if (access(filename, R_OK) != 0) {
+        perror(filename);
+        exit(1);
+    }
+}
+
+Options Options::parse(const int argc, char *const argv[]) {
+    string model;
+    string weights;
+    string means;
+    string dump;
+    string labels;
+
+    char c;
+
+    while ((c = getopt(argc, argv, "hm:w:n:l:d:")) != -1) {
+        switch (c) {
+            case 'h':
+                show_help(argv[0], 0);
+                break;
+
+            case 'w':
+                check_file(optarg);
+                weights = optarg;
+                break;
+
+            case 'n':
+                check_file(optarg);
+                model = optarg;
+                break;
+
+            case 'm':
+                check_file(optarg);
+                means = optarg;
+                break;
+
+            case 'l':
+                check_file(optarg);
+                labels = optarg;
+                break;
+
+            case 'd':
+                dump = optarg;
+                break;
+
+            case '?':
+                show_help(argv[0], 1);
+                break;
+
+            default:
+                cerr << "Unhandled option: " << c << endl;
+                abort();
+        }
+    }
+
+    if (weights.empty()) {
+        cerr << "Weights file is required!" << endl;
+        show_help(argv[0], 1);
+    }
+
+    if (model.empty()) {
+        cerr << "Model file is required!" << endl;
+        show_help(argv[0], 1);
+    }
+
+    for_each(argv + optind, argv + argc, check_file);
+
+    vector<string> inputs(argv + optind, argv + argc);
+    if (inputs.empty()) {
+        cerr << "No inputs specified" << endl;
+        show_help(argv[0], 1);
+    }
+
+    return Options(move(model), move(weights), move(means), move(labels), move(dump), move(inputs));
+}
+
+Options::Options(string &&model, string &&weights, string&& means, string&& labels, string&& dumpPath, vector<string> &&inputs) noexcept:
+        modelPath(move(model)),
+        weightsPath(move(weights)),
+        meansPath(means),
+        labelsPath(labels),
+        dumpPath(dumpPath),
+        inputPaths(move(inputs))
+{
+}
+
+const string& Options::model() const {
+    return modelPath;
+}
+
+const string& Options::weights() const {
+    return weightsPath;
+}
+
+const vector<string>& Options::inputs() const {
+    return inputPaths;
+}
+
+const string& Options::means() const
+{
+    return meansPath;
+}
+
+optional<vector<string>> Options::labels() const
+{
+    if (labelsPath.empty()) {
+        return nullopt;
+    } else {
+        return read_vector<string>(labelsPath);
+    }
+}
+
+std::optional<PNGDumper> Options::imageDumper() const
+{
+    if (dumpPath.empty()) {
+        return nullopt;
+    } else {
+        return move(PNGDumper(dumpPath));
+    }
+}

src/fmri/Options.hpp

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <optional>
+#include <string>
+#include <vector>
+
+#include "PNGDumper.hpp"
+
+namespace fmri {
+
+    using std::vector;
+    using std::string;
+
+    class Options {
+    public:
+        static Options parse(const int argc, char *const argv[]);
+
+        const string& model() const;
+        const string& weights() const;
+        const string& means() const;
+        std::optional<vector<string>> labels() const;
+        std::optional<fmri::PNGDumper> imageDumper() const;
+
+        const vector<string>& inputs() const;
+
+    private:
+        const string modelPath;
+        const string weightsPath;
+        const string meansPath;
+        const string labelsPath;
+        const string dumpPath;
+        const vector<string> inputPaths;
+
+        Options(string &&, string &&, string&&, string&&, string&&, vector<string> &&) noexcept;
+    };
+}

src/fmri/PNGDumper.cpp

@@ -0,0 +1,90 @@
+#include <cstring>
+
+#include <glog/logging.h>
+#include <sys/stat.h>
+#include <png++/png.hpp>
+
+#include "PNGDumper.hpp"
+
+using namespace fmri;
+using namespace std;
+
+static void ensureDir(const string& dir)
+{
+    struct stat s;
+    if (stat(dir.c_str(), &s) == 0) {
+        CHECK(S_ISDIR(s.st_mode)) << dir << " already exists and is not a directory." << endl;
+        return;
+    }
+
+    switch (errno) {
+        case ENOENT:
+            PCHECK(mkdir(dir.c_str(), 0777) == 0) << "Couldn't create directory";
+            return;
+
+        default:
+            perror("Unusable dump dir");
+            break;
+    }
+}
+
+PNGDumper::PNGDumper(string_view baseDir) :
+        baseDir_(baseDir)
+{
+    ensureDir(baseDir_);
+}
+
+void PNGDumper::dump(const LayerData &layerData)
+{
+    if (layerData.shape().size() == 4) {
+        // We have a series of images.
+        dumpImageSeries(layerData);
+    } else {
+        LOG(INFO) << "Unable to dump this type of layer to png.";
+    }
+}
+
+void PNGDumper::dumpImageSeries(const LayerData &layer)
+{
+    const auto& shape = layer.shape();
+    const auto images = shape[0], channels = shape[1], height = shape[2], width = shape[3];
+    const auto imagePixels = width * height;
+
+    // Buffer for storing the current image data.
+    vector<DType> buffer(imagePixels);
+
+    auto data = layer.data();
+
+    png::image<png::gray_pixel> image(width, height);
+
+    for (int i = 0; i < images; ++i) {
+        for (int j = 0; j < channels; ++j) {
+            memcpy(buffer.data(), data, imagePixels * sizeof(DType));
+
+            // advance the buffer;
+            data += imagePixels;
+
+            rescale(buffer.begin(), buffer.end(), 0.0, 255.0);
+
+            for (int y = 0; y < height; ++y) {
+                for (int x = 0; x < width; ++x) {
+                    image[y][x] = png::gray_pixel((int) buffer[x + y * width]);
+                }
+            }
+
+            image.write(getFilename(layer.name(), i, j));
+        }
+    }
+}
+
+string PNGDumper::getFilename(const string &layerName, int i, int j)
+{
+    stringstream nameBuilder;
+
+    nameBuilder << baseDir_
+                << "/" << layerName
+                << "-" << i
+                << "-" << j << ".png";
+
+    return nameBuilder.str();
+}

src/fmri/PNGDumper.hpp

@@ -0,0 +1,28 @@
+#pragma once
+
+#include <string>
+#include <string_view>
+
+#include "LayerData.hpp"
+#include "utils.hpp"
+
+namespace fmri
+{
+    using std::string;
+    using std::string_view;
+
+    class PNGDumper
+    {
+    public:
+        PNGDumper(string_view baseDir);
+
+        void dump(const LayerData& layerData);
+
+    private:
+        string baseDir_;
+
+        void dumpImageSeries(const LayerData &data);
+
+        string getFilename(const string &basic_string, int i, int j);
+    };
+}

src/fmri/PoolingLayerAnimation.cpp

@@ -0,0 +1,47 @@
+#include <glog/logging.h>
+#include <cmath>
+#include <caffe/util/math_functions.hpp>
+#include "PoolingLayerAnimation.hpp"
+#include "glutils.hpp"
+#include "MultiImageVisualisation.hpp"
+
+using namespace std;
+using namespace fmri;
+
+PoolingLayerAnimation::PoolingLayerAnimation(const LayerData &prevData, const LayerData &curData,
+                                             const std::vector<float> &prevPositions,
+                                             const std::vector<float> &curPositions) :
+        original(loadTextureForData(prevData)),
+        downSampled(loadTextureForData(curData)),
+        startingPositions(MultiImageVisualisation::getVertices(prevPositions)),
+        deltas(startingPositions.size()),
+        textureCoordinates(MultiImageVisualisation::getTexCoords(prevPositions.size() / 3))
+{
+    CHECK_EQ(prevPositions.size(), curPositions.size()) << "Layers should be same size. Caffe error?";
+    const auto downScaling = sqrt(
+            static_cast<float>(curData.shape()[2] * curData.shape()[3]) / (prevData.shape()[2] * prevData.shape()[3]));
+    const auto targetPositions = MultiImageVisualisation::getVertices(curPositions, downScaling);
+    caffe::caffe_sub(targetPositions.size(), targetPositions.data(), startingPositions.data(), deltas.data());
+
+    for (auto i = 0u; i < deltas.size(); i+=3) {
+        deltas[i] = LAYER_X_OFFSET;
+    }
+}
+
+void PoolingLayerAnimation::draw(float timeStep)
+{
+    vector<float> vertexBuffer(deltas);
+    caffe::caffe_scal(vertexBuffer.size(), timeStep, vertexBuffer.data());
+    caffe::caffe_add(startingPositions.size(), startingPositions.data(), vertexBuffer.data(), vertexBuffer.data());
+
+    drawImageTiles(vertexBuffer.size() / 3, vertexBuffer.data(), textureCoordinates.data(), original);
+}
+
+Texture PoolingLayerAnimation::loadTextureForData(const LayerData &data)
+{
+    CHECK_EQ(data.shape().size(), 4) << "Layer should be image-like";
+    CHECK_EQ(data.shape()[0], 1) << "Only single images supported";
+
+    auto channels = data.shape()[1], width = data.shape()[2], height = data.shape()[3];
+    return loadTexture(data.data(), width, height * channels, channels);
+}

src/fmri/PoolingLayerAnimation.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include "Animation.hpp"
+#include "LayerData.hpp"
+#include "Texture.hpp"
+
+namespace fmri
+{
+    class PoolingLayerAnimation : public Animation
+    {
+    public:
+        PoolingLayerAnimation(const LayerData &prevData, const LayerData &curData,
+                                      const std::vector<float> &prevPositions,
+                                      const std::vector<float> &curPositions);
+
+        void draw(float timeStep) override;
+
+    private:
+        Texture original;
+        Texture downSampled;
+        std::vector<float> startingPositions;
+        std::vector<float> deltas;
+        std::vector<float> textureCoordinates;
+
+        static Texture loadTextureForData(const LayerData& data);
+    };
+}

src/fmri/Range.hpp

@@ -0,0 +1,66 @@
+#pragma once
+
+namespace fmri
+{
+
+    /**
+     * Iterable that produces a specific range of integers.
+     *
+     * Useful to make an automatically typed for loop over an unknown integer type.
+     *
+     * @tparam T The integer type to use.
+     */
+    template<class T>
+    class Range
+    {
+    private:
+        T start_;
+        T end_;
+
+    public:
+        /**
+         * Construct a range from 0 to num, non inclusive.
+         *
+         * @param num
+         */
+        constexpr explicit Range(const T &num) : start_(0), end_(num) {};
+
+        /**
+         * Construct a range from start to end, non inclusive.
+         *
+         * @param start
+         * @param end
+         */
+        constexpr Range(const T &start, const T &end) : start_(start), end_(end) {};
+
+        class Iter
+        {
+        private:
+            T cur_;
+
+        public:
+            constexpr explicit Iter(const T &cur) : cur_(cur)
+            {};
+
+            typedef std::bidirectional_iterator_tag iterator_category;
+            typedef T value_type;
+            typedef typename std::make_signed<T>::type difference_type;
+            typedef T& reference;
+            typedef T* pointer;
+
+            constexpr bool operator!=(const Iter& o) { return o.cur_ != cur_; }
+            constexpr Iter&operator++() { ++cur_; return *this; }
+            constexpr Iter&operator--() { --cur_; return *this; }
+
+            constexpr const T &operator*() { return cur_; }
+        };
+
+        typedef Iter const_iterator;
+        typedef std::reverse_iterator<Iter> const_reverse_iterator;
+
+        constexpr const_iterator begin() const { return Iter(start_); }
+        constexpr const_iterator end() const { return Iter(end_); }
+        constexpr const_reverse_iterator rbegin() const { return std::make_reverse_iterator(end()); }
+        constexpr const_reverse_iterator rend() const { return std::make_reverse_iterator(begin()); }
+    };
+}

src/fmri/Simulator.cpp

@@ -0,0 +1,233 @@
+#include <cassert>
+#include <iostream>
+#include <optional>
+#include <vector>
+
+#include <caffe/caffe.hpp>
+#include <opencv2/core/core.hpp>
+#include <opencv2/highgui/highgui.hpp>
+#include <opencv2/imgproc/imgproc.hpp>
+
+#include "Simulator.hpp"
+#include "Range.hpp"
+
+using namespace caffe;
+using namespace std;
+using namespace fmri;
+
+struct Simulator::Impl
+{
+    caffe::Net<DType> net;
+    cv::Size input_geometry;
+    optional<cv::Mat> means;
+    unsigned int num_channels;
+    map<string, LayerInfo> layerInfo_;
+
+    Impl(const string& model_file, const string& weights_file, const string& means_file);
+
+    vector<cv::Mat> getWrappedInputLayer();
+    cv::Mat preprocess(cv::Mat original) const;
+    vector<LayerData> simulate(const string &input_file);
+    const map<string, LayerInfo>& layerInfo() const;
+
+    void computeLayerInfo();
+
+    void loadMeans(const string &means_file);
+
+    void ensureNoInPlaceLayers();
+};
+
+// Create simple forwarding functions.
+Simulator::Simulator(const string& model_file, const string& weights_file, const string& means_file) :
+	pImpl(new Impl(model_file, weights_file, means_file))
+{
+}
+
+vector<LayerData> Simulator::simulate(const string& image_file)
+{
+    return pImpl->simulate(image_file);
+}
+
+Simulator::Impl::Impl(const string& model_file, const string& weights_file, const string& means_file) :
+	net(model_file, TEST)
+{
+	net.CopyTrainedLayersFrom(weights_file);
+    ensureNoInPlaceLayers();
+
+	auto input_layer = net.input_blobs()[0];
+	input_geometry = cv::Size(input_layer->width(), input_layer->height());
+	num_channels = input_layer->channels();
+
+	input_layer->Reshape(1, num_channels,
+			input_geometry.height, input_geometry.width);
+	/* Forward dimension change to all layers. */
+	net.Reshape();
+
+    if (!means_file.empty()) {
+        loadMeans(means_file);
+    }
+
+    computeLayerInfo();
+}
+
+void Simulator::Impl::loadMeans(const string &means_file)
+{// Read in the means file
+    BlobProto proto;
+    ReadProtoFromBinaryFileOrDie(means_file, &proto);
+
+    Blob<DType> mean_blob;
+    mean_blob.FromProto(proto);
+
+    CHECK_EQ(mean_blob.channels(), num_channels) << "Number of channels should match!" << endl;
+
+    vector<cv::Mat> channels;
+    float *data = mean_blob.mutable_cpu_data();
+    for (auto i : Range(num_channels)) {
+        (void)i;// Suppress unused warning
+        channels.emplace_back(mean_blob.height(), mean_blob.width(), CV_32FC1, data);
+        data += mean_blob.height() * mean_blob.width();
+    }
+
+    cv::Mat mean;
+    merge(channels, mean);
+
+    this->means = cv::Mat(input_geometry, mean.type(), cv::mean(mean));
+}
+
+vector<LayerData> Simulator::Impl::simulate(const string& image_file)
+{
+	cv::Mat im = cv::imread(image_file, -1);
+
+    assert(!im.empty());
+
+    auto input = preprocess(im);
+    auto channels = getWrappedInputLayer();
+
+    cv::split(input, channels);
+
+    net.Forward();
+
+	vector<LayerData> result;
+
+	const auto& names = net.layer_names();
+	const auto& results = net.top_vecs();
+
+	for (auto i : Range(names.size())) {
+		CHECK_EQ(results[i].size(), 1) << "Multiple outputs per layer are not supported!" << endl;
+		const auto blob = results[i][0];
+
+		result.emplace_back(names[i], blob->shape(), blob->cpu_data());
+	}
+
+    return result;
+}
+
+vector<cv::Mat> Simulator::Impl::getWrappedInputLayer()
+{
+    vector<cv::Mat> channels;
+    auto input_layer = net.input_blobs()[0];
+
+    const int width = input_geometry.width;
+    const int height = input_geometry.height;
+
+    DType* input_data = input_layer->mutable_cpu_data();
+    for (auto i : Range(num_channels)) {
+        (void)i;// Suppress unused warning
+        channels.emplace_back(height, width, CV_32FC1, input_data);
+        input_data += width * height;
+    }
+
+    return channels;
+}
+
+static cv::Mat fix_channels(const int num_channels, cv::Mat original) {
+    cv::Mat converted;
+
+    if (num_channels == 1 && original.channels() == 3) {
+        cv::cvtColor(original, converted, cv::COLOR_BGR2GRAY);
+    } else if (num_channels == 1 && original.channels() == 4) {
+        cv::cvtColor(original, converted, cv::COLOR_BGRA2GRAY);
+    } else if (num_channels == 3 && original.channels() == 1) {
+        cv::cvtColor(original, converted, cv::COLOR_GRAY2BGR);
+    } else if (num_channels == 3 && original.channels() == 4) {
+        cv::cvtColor(original, converted, cv::COLOR_BGRA2BGR);
+    } else {
+		CHECK(num_channels == original.channels()) << "Cannot convert between channel types. ";
+        return original;
+    }
+
+    return converted;
+}
+
+static cv::Mat resize(const cv::Size& targetSize, cv::Mat original)
+{
+    if (targetSize != original.size()) {
+        cv::Mat resized;
+        cv::resize(original, resized, targetSize);
+
+        return resized;
+    }
+
+    return original;
+}
+
+cv::Mat Simulator::Impl::preprocess(cv::Mat original) const
+{
+    auto converted = fix_channels(num_channels, std::move(original));
+
+    auto resized = resize(input_geometry, converted);
+
+    cv::Mat sample_float;
+    resized.convertTo(sample_float, num_channels == 3 ? CV_32FC3 : CV_32FC1);
+
+    if (!means) {
+        return sample_float;
+    }
+
+    cv::Mat normalized;
+    cv::subtract(sample_float, *means, normalized);
+
+    return normalized;
+
+}
+
+const map<string, LayerInfo> &Simulator::Impl::layerInfo() const
+{
+    return layerInfo_;
+}
+
+void Simulator::Impl::computeLayerInfo()
+{
+    const auto& names = net.layer_names();
+    const auto& layers = net.layers();
+
+    CHECK_EQ(names.size(), layers.size()) << "Size mismatch";
+
+    for (auto i : Range(names.size())) {
+        auto& layer = layers[i];
+        LayerInfo layerInfo(names[i], layer->type(), layer->blobs());
+        CHECK_NE(layerInfo.type(), LayerInfo::Type::Split) << "Split layers are not supported!";
+        layerInfo_.emplace(names[i], std::move(layerInfo));
+    }
+}
+
+void Simulator::Impl::ensureNoInPlaceLayers()
+{
+    auto blobList = net.top_vecs();
+    typeof(blobList) uniqueVecs;
+    unique_copy(blobList.begin(), blobList.end(), back_inserter(uniqueVecs));
+
+    LOG_IF(ERROR, blobList.size() != uniqueVecs.size())
+        << "Network file contains in-place layers, layer-state will not be accurate\n"
+        << "If accurate results are desired, see the deinplace script in tools." << endl;
+}
+
+Simulator::~Simulator()
+{
+    // Empty but defined constructor.
+}
+
+const map<string, LayerInfo> & Simulator::layerInfo() const
+{
+    return pImpl->layerInfo();
+}

src/fmri/Simulator.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <string>
+#include <memory>
+#include <vector>
+#include "LayerData.hpp"
+#include "LayerInfo.hpp"
+
+namespace fmri {
+    using std::string;
+    using std::vector;
+
+    class Simulator {
+    public:
+        Simulator(const string &model_file, const string &weights_file, const string &means_file = "");
+        ~Simulator();
+
+        vector<LayerData> simulate(const string &input_file);
+		const std::map<std::string, LayerInfo>& layerInfo() const;
+
+    private:
+		struct Impl;
+		std::unique_ptr<Impl> pImpl;
+    };
+}

src/fmri/Texture.cpp

@@ -0,0 +1,51 @@
+#include <algorithm>
+#include "Texture.hpp"
+
+using namespace fmri;
+
+Texture::Texture() noexcept
+{
+    glGenTextures(1, &id);
+}
+
+Texture::Texture(GLuint id) noexcept : id(id)
+{
+}
+
+Texture::~Texture()
+{
+    if (id != 0) {
+        glDeleteTextures(1, &id);
+    }
+}
+
+Texture &Texture::operator=(Texture && other) noexcept
+{
+    std::swap(id, other.id);
+    return *this;
+}
+
+Texture::Texture(Texture && other) noexcept
+{
+    std::swap(id, other.id);
+}
+
+void Texture::bind(GLenum target) const
+{
+    glBindTexture(target, id);
+}
+
+void Texture::configure(GLenum target) const
+{
+    bind(target);
+    const float color[] = {1, 0, 1}; // Background color for textures.
+
+    // Set up (lack of) repetition
+    glTexParameteri(target, GL_TEXTURE_WRAP_R, GL_CLAMP_TO_BORDER);
+    glTexParameteri(target, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_BORDER);
+    glTexParameterfv(target, GL_TEXTURE_BORDER_COLOR, color);
+
+    // Set up texture scaling
+    glTexParameteri(target, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST); // Use mipmapping for scaling down
+    glTexParameteri(target, GL_TEXTURE_MAG_FILTER, GL_NEAREST); // Use nearest pixel when scaling up.
+}

src/fmri/Texture.hpp

@@ -0,0 +1,45 @@
+#pragma once
+
+#include <GL/gl.h>
+
+namespace fmri
+{
+    /**
+     * Simple owning Texture class.
+     *
+     * Encapsulates an OpenGL texture, and enables RAII for it. Copying
+     * is disallowed for this reason.
+     */
+    class Texture
+    {
+    public:
+        /**
+         * Allocate a new texture
+         */
+        Texture() noexcept;
+        Texture(Texture &&) noexcept;
+        Texture(const Texture &) = delete;
+
+        /**
+         * Own an existing texture
+         * @param id original texture ID.
+         */
+        explicit Texture(GLuint id) noexcept;
+
+        ~Texture();
+
+        Texture &operator=(Texture &&) noexcept;
+        Texture &operator=(const Texture &) = delete;
+
+        /**
+         * Bind the owned texture to the given spot.
+         * @param target valid target for glBindTexture.
+         */
+        void bind(GLenum target) const;
+        void configure(GLenum target) const;
+
+    private:
+        GLuint id;
+
+    };
+}

src/fmri/camera.cpp

@@ -0,0 +1,134 @@
+#include <GL/freeglut.h>
+#include <cmath>
+#include <sstream>
+#include <iostream>
+#include "camera.hpp"
+#include "utils.hpp"
+
+using namespace fmri;
+using namespace std;
+
+static Camera& camera = Camera::instance();
+
+static void handleMouseMove(int x, int y)
+{
+    const float width = glutGet(GLUT_WINDOW_WIDTH) / 2;
+    const float height = glutGet(GLUT_WINDOW_HEIGHT) / 2;
+
+    camera.angle[0] = (x - width) / width * 180;
+    camera.angle[1] = (y - height) / height * 90;
+
+    glutPostRedisplay();
+}
+
+static float getFPS()
+{
+    static int frames = 0;
+    static float fps = 0;
+    static auto timeBase = glutGet(GLUT_ELAPSED_TIME);
+
+    ++frames;
+    const auto time = glutGet(GLUT_ELAPSED_TIME);
+    if (time - timeBase > 2000) {
+        fps = frames * 1000.0f / (time - timeBase);
+        frames = 0;
+        timeBase = time;
+    }
+
+
+    return fps;
+}
+
+static void move(unsigned char key)
+{
+    float speed = 0.5;
+    float dir[3];
+
+    const auto yaw = deg2rad(camera.angle[0]);
+    const auto pitch = deg2rad(camera.angle[1]);
+
+    if (key == 'w' || key == 's') {
+        dir[0] = sin(yaw) * cos(pitch);
+        dir[1] = -sin(pitch);
+        dir[2] = -cos(yaw) * cos(pitch);
+    } else {
+        dir[0] = -cos(yaw);
+        dir[1] = 0;
+        dir[2] = -sin(yaw);
+    }
+
+    if (key == 's' || key == 'd') {
+        speed *= -1;
+    }
+
+    for (auto i = 0; i < 3; ++i) {
+        camera.pos[i] += speed * dir[i];
+    }
+}
+
+static void handleKeys(unsigned char key, int, int)
+{
+    switch (key) {
+        case 'w':
+        case 'a':
+        case 's':
+        case 'd':
+            move(key);
+            break;
+
+        case 'q':
+            // Utility quit function.
+            glutLeaveMainLoop();
+            break;
+
+        case 'h':
+            camera.reset();
+            break;
+
+        default:
+            // Do nothing.
+            break;
+    }
+
+    glutPostRedisplay();
+}
+
+std::string Camera::infoLine()
+{
+    stringstream buffer;
+    buffer << "Pos(x,y,z) = (" << pos[0] << ", " << pos[1] << ", " << pos[2] << ")\n";
+    buffer << "Angle(p,y) = (" << angle[0] << ", " << angle[1] << ")\n";
+    buffer << "FPS = " << getFPS() << "\n";
+    return buffer.str();
+}
+
+void Camera::reset()
+{
+    pos[0] = 0;
+    pos[1] = 0;
+    pos[2] = 3;
+
+    angle[0] = 0;
+    angle[1] = 0;
+}
+
+void Camera::configureRenderingContext()
+{
+    glLoadIdentity();
+    glRotatef(angle[1], 1, 0, 0);
+    glRotatef(angle[0], 0, 1, 0);
+    glTranslatef(-pos[0], -pos[1], -pos[2]);
+}
+
+Camera &Camera::instance()
+{
+    static Camera camera;
+    return camera;
+}
+
+void Camera::registerControls()
+{
+    reset();
+    glutPassiveMotionFunc(handleMouseMove);
+    glutKeyboardFunc(handleKeys);
+}

src/fmri/camera.hpp

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <string>
+
+namespace fmri
+{
+    struct Camera {
+        float pos[3];
+        float angle[2];
+
+        void reset();
+        void configureRenderingContext();
+        void registerControls();
+        std::string infoLine();
+
+        static Camera& instance();
+
+    private:
+        Camera() noexcept = default;
+    };
+}

src/fmri/glutils.cpp

@@ -0,0 +1,154 @@
+#include <cstdint>
+#include <memory>
+#include <vector>
+#include <cstring>
+#include <glog/logging.h>
+#include <chrono>
+#include <thread>
+#include "glutils.hpp"
+#include "Range.hpp"
+
+using namespace fmri;
+using namespace std;
+
+static void handleGLError(GLenum error) {
+    switch (error) {
+        case GL_NO_ERROR:
+            return;
+
+        default:
+            cerr << "OpenGL error: " << (const char*) gluGetString(error) << endl;
+    }
+}
+
+static void rescaleSubImages(vector<float>& textureBuffer, int subImages) {
+    auto cur = textureBuffer.begin();
+    const auto increment = textureBuffer.size() / subImages;
+
+    while (cur != textureBuffer.end()) {
+        rescale(cur, cur + increment, 0, 1);
+        advance(cur, increment);
+    }
+}
+
+fmri::Texture fmri::loadTexture(DType const *data, int width, int height, int subImages)
+{
+    // Load and scale texture
+    vector<float> textureBuffer(data, data + (width * height));
+    rescaleSubImages(textureBuffer, subImages);
+
+    Texture texture;
+    texture.configure(GL_TEXTURE_2D);
+    gluBuild2DMipmaps(GL_TEXTURE_2D, GL_LUMINANCE, width, height, GL_LUMINANCE, GL_FLOAT, textureBuffer.data());
+
+    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, 10000.0f);
+
+    // Get Back to the Modelview
+    glMatrixMode(GL_MODELVIEW);
+}
+
+void fmri::renderText(std::string_view text, int x, int y)
+{
+    constexpr auto font = GLUT_BITMAP_HELVETICA_10;
+    glRasterPos2i(x, y);
+    for (char c : text) {
+        if (c == '\n') {
+            y += 12;
+            glRasterPos2i(x, y);
+        } else {
+            glutBitmapCharacter(font, c);
+        }
+    }
+}
+
+void fmri::checkGLErrors()
+{
+    while (auto error = glGetError()) {
+        handleGLError(error);
+    }
+}
+
+void fmri::throttleIdleFunc()
+{
+    using namespace std::chrono;
+
+    constexpr duration<double, ratio<1, 60>> refreshRate(1);
+
+    static auto lastCalled = steady_clock::now();
+
+    const auto now = steady_clock::now();
+
+    const auto diff = now - lastCalled;
+
+    if (diff < refreshRate) {
+        this_thread::sleep_for(refreshRate - diff);
+    }
+
+    lastCalled = now;
+}
+
+void fmri::restorePerspectiveProjection() {
+
+    glMatrixMode(GL_PROJECTION);
+    // restore previous projection matrix
+    glPopMatrix();
+
+    // get back to modelview mode
+    glMatrixMode(GL_MODELVIEW);
+}
+
+void fmri::setOrthographicProjection() {
+
+    // switch to projection mode
+    glMatrixMode(GL_PROJECTION);
+
+    // save previous matrix which contains the
+    //settings for the perspective projection
+    glPushMatrix();
+
+    // reset matrix
+    glLoadIdentity();
+
+    // set a 2D orthographic projection
+    gluOrtho2D(0, glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT), 0);
+
+    // switch back to modelview mode
+    glMatrixMode(GL_MODELVIEW);
+}
+
+void fmri::drawImageTiles(int n, const float *vertexBuffer, const float *textureCoords, const Texture &texture)
+{
+    glEnableClientState(GL_TEXTURE_COORD_ARRAY);
+    glEnableClientState(GL_VERTEX_ARRAY);
+    glEnable(GL_TEXTURE_2D);
+    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
+    texture.bind(GL_TEXTURE_2D);
+    glTexCoordPointer(2, GL_FLOAT, 0, textureCoords);
+    glVertexPointer(3, GL_FLOAT, 0, vertexBuffer);
+    glDrawArrays(GL_QUADS, 0, n);
+    glDisable(GL_TEXTURE_2D);
+    glDisableClientState(GL_VERTEX_ARRAY);
+    glDisableClientState(GL_TEXTURE_COORD_ARRAY);
+}

src/fmri/glutils.hpp

@@ -0,0 +1,64 @@
+#pragma once
+
+#include "LayerData.hpp"
+#include "utils.hpp"
+#include "Texture.hpp"
+#include <GL/glut.h>
+#include <string_view>
+
+namespace fmri {
+    /**
+     * Create a texture from an array of data.
+     *
+     * @param data
+     * @param width
+     * @param height
+     * @param subImages Number of subimages in the original image. Sub images are rescaled individually to preserve contrast. Optional, default 1.
+     * @return A texture reference.
+     */
+    fmri::Texture loadTexture(DType const *data, int width, int height, int subImages = 1);
+
+    /**
+     * 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);
+
+    /**
+     * Draw a bitmap string at the current location.
+     *
+     * @param text The text to draw.
+     */
+    void renderText(std::string_view text, int x = 0, int y = 0);
+
+    /**
+     * Check if there are OpenGL errors and report them.
+     */
+    void checkGLErrors();
+
+    /**
+     * Slow down until the idle func is being called a reasonable amount of times.
+     */
+    void throttleIdleFunc();
+
+    void setOrthographicProjection();
+
+    void restorePerspectiveProjection();
+
+    /**
+     * Draw a series of textured tiles to the screen.
+     *
+     * This function ends up drawing GL_QUADS.
+     *
+     * @param n Number of vertices
+     * @param vertexBuffer
+     * @param textureCoords
+     * @param texture
+     */
+    void drawImageTiles(int n, const float* vertexBuffer, const float* textureCoords, const Texture& texture);
+}

src/fmri/main.cpp

@@ -0,0 +1,216 @@
+#include <algorithm>
+#include <iostream>
+#include <glog/logging.h>
+#include <GL/glut.h>
+#include <map>
+
+#include "LayerData.hpp"
+#include "Options.hpp"
+#include "Simulator.hpp"
+#include "glutils.hpp"
+#include "camera.hpp"
+#include "LayerVisualisation.hpp"
+#include "Range.hpp"
+#include "visualisations.hpp"
+
+using namespace std;
+using namespace fmri;
+
+struct
+{
+    optional<vector<string>> labels;
+    map<string, LayerInfo> layerInfo;
+    vector<vector<LayerData>> data;
+    vector<vector<LayerData>>::iterator currentData;
+    vector<unique_ptr<LayerVisualisation>> layerVisualisations;
+    vector<unique_ptr<Animation>> animations;
+    float animationStep = 0;
+} rendererData;
+
+static void loadSimulationData(const Options &options)
+{
+    vector<vector<LayerData>> &results = rendererData.data;
+    results.clear();
+
+    auto dumper = options.imageDumper();
+    Simulator simulator(options.model(), options.weights(), options.means());
+    rendererData.layerInfo = simulator.layerInfo();
+
+    for (const auto &image : options.inputs()) {
+        results.emplace_back(simulator.simulate(image));
+    }
+
+    CHECK_GT(results.size(), 0) << "Should have some results" << endl;
+
+    if (dumper) {
+        for (auto &layer : *results.begin()) {
+            dumper->dump(layer);
+        }
+    }
+
+    const auto optLabels = options.labels();
+
+    if (optLabels) {
+        auto& labels = *optLabels;
+        for (const auto& result : results) {
+            auto& last = *result.rbegin();
+            auto bestIndex = std::distance(last.data(), max_element(last.data(), last.data() + last.numEntries()));
+            LOG(INFO) << "Got answer: " << labels[bestIndex] << endl;
+        }
+    }
+}
+
+static void renderLayerName(const LayerData &data);
+
+static void renderDebugInfo()
+{
+    glLoadIdentity();
+    setOrthographicProjection();
+    glColor3f(1, 1, 0);
+    renderText(Camera::instance().infoLine(), 2, 10);
+    restorePerspectiveProjection();
+}
+
+static void render()
+{
+    // Clear Color and Depth Buffers
+    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+    auto& camera = Camera::instance();
+
+    camera.configureRenderingContext();
+
+    const auto& dataSet = *rendererData.currentData;
+
+    glPushMatrix();
+    glTranslatef(5 * dataSet.size(), 0, 0);
+
+    for (auto i : Range(dataSet.size())) {
+        glPushMatrix();
+        renderLayerName(dataSet[i]);
+        rendererData.layerVisualisations[i]->render();
+        if (i < rendererData.animations.size() && rendererData.animations[i]) {
+            rendererData.animations[i]->draw(rendererData.animationStep);
+        }
+
+        glPopMatrix();
+        glTranslatef(LAYER_X_OFFSET, 0, 0);
+    }
+
+
+    glPopMatrix();
+
+    renderDebugInfo();
+
+    glutSwapBuffers();
+}
+
+static void renderLayerName(const LayerData &data)
+{
+    // Draw the name of the layer for reference.
+    glColor3f(0.5, 0.5, 0.5);
+    renderText(data.name());
+
+    glTranslatef(0, 0, -10);
+}
+
+static void updateVisualisers()
+{
+    rendererData.layerVisualisations.clear();
+    rendererData.animations.clear();
+    LayerData* prevState = nullptr;
+    LayerVisualisation* prevVisualisation = nullptr;
+
+    for (LayerData &layer : *rendererData.currentData) {
+        LayerVisualisation* visualisation = getVisualisationForLayer(layer, rendererData.layerInfo.at(layer.name()));
+        if (prevState && prevVisualisation && visualisation) {
+            auto interaction = getActivityAnimation(*prevState, layer, rendererData.layerInfo.at(layer.name()), prevVisualisation->nodePositions(), visualisation->nodePositions());
+            rendererData.animations.emplace_back(interaction);
+        }
+
+        rendererData.layerVisualisations.emplace_back(visualisation);
+
+        prevVisualisation = visualisation;
+        prevState = &layer;
+    }
+
+    glutPostRedisplay();
+}
+
+static void specialKeyFunc(int key, int, int)
+{
+    switch (key) {
+        case GLUT_KEY_LEFT:
+            if (rendererData.currentData == rendererData.data.begin()) {
+                rendererData.currentData = rendererData.data.end();
+            }
+            --rendererData.currentData;
+            updateVisualisers();
+            break;
+
+        case GLUT_KEY_RIGHT:
+            ++rendererData.currentData;
+            if (rendererData.currentData == rendererData.data.end()) {
+                rendererData.currentData = rendererData.data.begin();
+            }
+            updateVisualisers();
+            break;
+
+        default:
+            LOG(INFO) << "Received keystroke " << key;
+    }
+}
+
+static void idleFunc()
+{
+    checkGLErrors();
+    glutPostRedisplay();
+    throttleIdleFunc();
+    rendererData.animationStep = (1 - cos(M_PI * getAnimationStep(std::chrono::seconds(5)))) / 2;
+}
+
+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();
+    loadSimulationData(options);
+
+    // Register callbacks
+    glutDisplayFunc(render);
+    glutIdleFunc(idleFunc);
+    glutReshapeFunc(changeWindowSize);
+    glutSpecialFunc(specialKeyFunc);
+
+    Camera::instance().registerControls();
+
+    rendererData.currentData = rendererData.data.begin();
+    updateVisualisers();
+
+    // Enable depth test to fix objects behind you
+    glEnable(GL_DEPTH_TEST);
+
+    // Nicer rendering
+    glEnable(GL_POINT_SMOOTH);
+    glEnable(GL_LINE_SMOOTH);
+    glEnable(GL_POLYGON_SMOOTH);
+    glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);
+    glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);
+    glHint(GL_POLYGON_SMOOTH_HINT, GL_NICEST);
+    glEnable(GL_BLEND);
+    glBlendFunc(GL_DST_ALPHA,GL_ONE_MINUS_DST_ALPHA);
+
+    // Start visualisation
+    glutMainLoop();
+
+    google::ShutdownGoogleLogging();
+
+    return 0;
+}

src/fmri/utils.cpp

@@ -0,0 +1,16 @@
+#include "utils.hpp"
+
+const float fmri::LAYER_X_OFFSET = -10;
+
+std::default_random_engine &fmri::rng()
+{
+    static std::default_random_engine rng;
+    static std::default_random_engine::result_type seed = 0;
+
+    if (seed == 0) {
+        std::random_device dev;
+        rng.seed(seed = dev());
+    }
+
+    return rng;
+}

src/fmri/utils.hpp

@@ -0,0 +1,232 @@
+#pragma once
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <fstream>
+#include <iterator>
+#include <random>
+#include <string>
+#include <utility>
+#include <vector>
+#include <ratio>
+#include <chrono>
+
+namespace fmri
+{
+	typedef float DType;
+
+    /**
+     * The distance between layers in the visualisation.
+     */
+	extern const float LAYER_X_OFFSET;
+
+    /**
+     * Identity function that simply returns whatever is put in.
+     *
+     * @tparam T The type of the function
+     * @param t The value to return.
+     * @return The original value.
+     */
+    template<class T>
+    inline T identity(T t) {
+        return t;
+    }
+
+    /**
+     * Read a file into a vector of given type.
+     * @tparam T The type to read.
+     * @param filename the file to read
+     * @return A vector of type T.
+     */
+    template<class T>
+    inline std::vector <T> read_vector(const std::string& filename)
+    {
+        std::ifstream input(filename);
+        assert(input.good());
+
+        std::vector<T> res;
+        std::transform(std::istream_iterator<T>(input),
+                       std::istream_iterator<T>(),
+                       identity<T>, std::back_inserter(res));
+
+        return res;
+    }
+
+    /**
+     * String specialisation of read_vector.
+     *
+     * @param filename The filename to load.
+     * @return A vector of the lines in the source file.
+     */
+    template<>
+    inline std::vector<std::string> read_vector<std::string>(const std::string& filename)
+    {
+        std::ifstream input(filename);
+        assert(input.good());
+
+        std::string v;
+        std::vector<std::string> res;
+
+        while (getline(input, v)) {
+            res.push_back(v);
+        }
+
+        return res;
+    }
+
+    /**
+     * Create a vector of pairs.
+     *
+     * @tparam T The first type
+     * @tparam U The second type
+     * @param a First vector
+     * @param b Second vector
+     * @return A vector of pair<U, V>
+     */
+    template<class T, class U>
+    std::vector<std::pair<T, U>> combine(const std::vector<T>& a, const std::vector<U>& b)
+    {
+        std::vector<std::pair<T, U>> res;
+        std::transform(a.begin(), a.end(), b.begin(),
+                       std::back_inserter(res),
+                       [] (const T& a, const U& b) -> auto { return std::make_pair(a, b); });
+
+        return res;
+    }
+
+    /**
+     * @brief Scales a range of values into a fixed range of values.
+     *
+     * This method traverses the range twice, once to determine maximum
+     * and minimum, and once again to modify all the values.
+     *
+     * @tparam It Iterator type. Will be used to determine value type.
+     * @param first The start of the range to scale
+     * @param last The end of the range to scale
+     * @param minimum The desired minimum of the range.
+     * @param maximum The desired maximum of the range.
+     */
+    template<class It>
+    void rescale(const It first, const It last,
+                 typename std::iterator_traits<It>::value_type minimum,
+                 typename std::iterator_traits<It>::value_type maximum)
+    {
+        const auto[minElem, maxElem] = std::minmax_element(first, last);
+        const auto rangeWidth = maximum - minimum;
+        const auto valWidth = *maxElem - *minElem;
+
+        if (valWidth == 0) {
+            // Just fill the range with the minimum value, since
+            std::fill(first, last, minimum);
+        } else {
+            const auto scaling = rangeWidth / valWidth;
+
+            const auto minVal = *minElem;
+
+            std::for_each(first, last, [=](typename std::iterator_traits<It>::reference v) {
+                v = std::clamp((v - minVal) * scaling + minimum, minimum, maximum);
+            });
+        }
+    }
+
+    template<class T>
+    constexpr inline T deg2rad(T val) {
+        return val / 180 * M_PI;
+    }
+
+    /**
+     * Compute the ideal number of columns for dividing up a range into a rectangle.
+     *
+     * @tparam T
+     * @param elems number of elements to space out
+     * @return the ideal number of columns
+     */
+    template<class T>
+    inline T numCols(const T elems) {
+        auto cols = static_cast<T>(ceil(sqrt(elems)));
+
+        while (elems % cols) {
+            ++cols;// TODO: this should probably be done analytically
+        }
+
+        return cols;
+    }
+
+    /**
+     * Get a globally initialized random number generator.
+     *
+     * This RNG should always be used as a reference, to make sure the state actually updates.
+     *
+     * @return A reference to the global RNG.
+     */
+    std::default_random_engine& rng();
+
+    /**
+     * Get the current animation offset for a particular animation.
+     *
+     * @tparam Duration Duration type of length. Should be a specialisation of std::chrono::duration
+     * @param length The length of the animation.
+     * @return
+     */
+    template<class Duration>
+    float getAnimationStep(const Duration &length) {
+        using namespace std::chrono;
+
+        static auto startingPoint = steady_clock::now();
+        const auto modified_length = duration_cast<steady_clock::duration>(length);
+
+        auto step = (steady_clock::now() - startingPoint) % modified_length;
+
+        return static_cast<float>(step.count()) / static_cast<float>(modified_length.count());
+    }
+
+    /**
+     * Perform an argsort partitioning on the first n elements.
+     *
+     * @tparam Iter
+     * @tparam Compare
+     * @param first First element
+     * @param middle Sorting limit
+     * @param last Past end iterator for range
+     * @param compare Comparison function to use
+     * @return A vector of the indices before the partitioning cut-off.
+     */
+    template<class Iter, class Compare>
+    std::vector<std::size_t> arg_nth_element(Iter first, Iter middle, Iter last, Compare compare)
+    {
+        using namespace std;
+
+        const auto n = static_cast<size_t>(distance(first, middle));
+        const auto total = static_cast<size_t>(distance(first, last));
+
+        vector<size_t> indices(total);
+        iota(indices.begin(), indices.end(), 0u);
+
+        nth_element(indices.begin(), indices.begin() + n, indices.end(), [=](size_t a, size_t b) {
+            return compare(*(first + a), *(first + b));
+        });
+
+        indices.resize(n);
+
+        return indices;
+    }
+
+    /**
+     * Fix non-normal floating point values in a range.
+     *
+     * @tparam It
+     * @param first Start of range iterator
+     * @param last Past the end of range iterator
+     * @param normalValue Value to assign to non-normal values. Default 1.
+     */
+    template<class It>
+    inline void normalize(It first, It last, typename std::iterator_traits<It>::value_type normalValue = 1)
+    {
+        for (; first != last; ++first) {
+            if (!std::isnormal(*first)) {
+                *first = normalValue;
+            }
+        }
+    }
+}

src/fmri/visualisations.cpp

@@ -0,0 +1,245 @@
+#include <algorithm>
+#include <numeric>
+#include <caffe/util/math_functions.hpp>
+#include "visualisations.hpp"
+#include "DummyLayerVisualisation.hpp"
+#include "MultiImageVisualisation.hpp"
+#include "FlatLayerVisualisation.hpp"
+#include "Range.hpp"
+#include "ActivityAnimation.hpp"
+#include "InputLayerVisualisation.hpp"
+#include "PoolingLayerAnimation.hpp"
+#include "ImageInteractionAnimation.hpp"
+
+using namespace fmri;
+using namespace std;
+
+// Maximum number of interactions shown
+static constexpr size_t INTERACTION_LIMIT = 10000;
+
+typedef vector<pair<float, pair<size_t, size_t>>> EntryList;
+
+/**
+ * Normalizer for node positions.
+ *
+ * Since not every neuron in a layer may get a node in the visualisation,
+ * this function maps those neurons back to a node number that does.
+ *
+ * Usage: node / getNodeNormalizer(layer).
+ *
+ * @param layer Layer to compute normalization for
+ * @return Number to divide node numbers by.
+ */
+static inline int getNodeNormalizer(const LayerData& layer) {
+    const auto& shape = layer.shape();
+    switch(shape.size()) {
+        case 2:
+            return 1;
+
+        case 4:
+            return shape[2] * shape[3];
+
+        default:
+            CHECK(false) << "Unsupported shape " << shape.size() << endl;
+            exit(EINVAL);
+    }
+}
+
+/**
+ * Deduplicate interaction entries.
+ *
+ * For duplicate interactions, the interaction strengths are summed.
+ *
+ * @param entries
+ * @return the deduplicated entries.
+ */
+static EntryList deduplicate(const EntryList& entries)
+{
+    map<pair<size_t, size_t>, float> combiner;
+    for (auto entry : entries) {
+        combiner[entry.second] += entry.first;
+    }
+
+    EntryList result;
+    transform(combiner.begin(), combiner.end(), back_inserter(result), [](const auto& item) {
+        return make_pair(item.second, item.first);
+    });
+
+    return result;
+}
+
+fmri::LayerVisualisation *fmri::getVisualisationForLayer(const fmri::LayerData &data, const fmri::LayerInfo &info)
+{
+    switch (info.type()) {
+        case LayerInfo::Type::Input:
+            if (data.shape().size() == 4) {
+                return new InputLayerVisualisation(data);
+            } else {
+                return new FlatLayerVisualisation(data, FlatLayerVisualisation::Ordering::SQUARE);
+            }
+
+        default:
+            switch (data.shape().size()) {
+                case 2:
+                    return new FlatLayerVisualisation(data, FlatLayerVisualisation::Ordering::SQUARE);
+
+                case 4:
+                    return new MultiImageVisualisation(data);
+
+                default:
+                    return new DummyLayerVisualisation();
+            }
+    }
+}
+
+static Animation *getFullyConnectedAnimation(const fmri::LayerData &prevState, const fmri::LayerInfo &layer,
+                                             const vector<float> &prevPositions, const vector<float> &curPositions)
+{
+    LOG(INFO) << "Computing top interactions for " << layer.name() << endl;
+
+    auto data = prevState.data();
+
+    CHECK_GE(layer.parameters().size(), 1) << "Layer should have correct parameters";
+
+    const auto shape = layer.parameters()[0]->shape();
+    auto weights = layer.parameters()[0]->cpu_data();
+    const auto numEntries = accumulate(shape.begin(), shape.end(), static_cast<size_t>(1), multiplies<void>());
+
+    vector<float> interactions(numEntries);
+    const auto stepSize = shape[0];
+
+    for (auto i : Range(numEntries / stepSize)) {
+        caffe::caffe_mul(shape[0], &weights[i * stepSize], data, &interactions[i * stepSize]);
+    }
+
+    const auto desiredSize = min(INTERACTION_LIMIT, numEntries);
+    auto idx = arg_nth_element(interactions.begin(), interactions.begin() + desiredSize, interactions.end(), [](auto a, auto b) {
+        return abs(a) > abs(b);
+    });
+
+    EntryList result;
+    result.reserve(desiredSize);
+    const auto normalizer = getNodeNormalizer(prevState);
+    for (auto i : idx) {
+        result.emplace_back(interactions[i], make_pair(i / shape[0] / normalizer, i % shape[0]));
+    }
+
+    return new ActivityAnimation(result, prevPositions.data(), curPositions.data());
+}
+
+static Animation *getDropOutAnimation(const fmri::LayerData &prevState,
+                                      const fmri::LayerData &curState,
+                                      const vector<float> &prevPositions,
+                                      const vector<float> &curPositions) {
+    const auto sourceNormalize = getNodeNormalizer(prevState);
+    const auto sinkNormalize = getNodeNormalizer(curState);
+
+    auto data = curState.data();
+    EntryList results;
+    results.reserve(curState.numEntries());
+    for (auto i : Range(curState.numEntries())) {
+        if (data[i] != 0) {
+            results.emplace_back(data[i], make_pair(i / sourceNormalize, i / sinkNormalize));
+        }
+    }
+
+    results = deduplicate(results);
+
+    return new ActivityAnimation(results, prevPositions.data(), curPositions.data());
+}
+
+static Animation *getReLUAnimation(const fmri::LayerData &prevState,
+                                   const fmri::LayerData &curState,
+                                   const vector<float> &prevPositions,
+                                   const vector<float> &curPositions) {
+    CHECK_EQ(curState.numEntries(), prevState.numEntries()) << "Layers should be of same size!";
+
+    std::vector<float> changes(prevState.numEntries());
+    caffe::caffe_sub(prevState.numEntries(), curState.data(), prevState.data(), changes.data());
+
+    if (curState.shape().size() == 2) {
+        EntryList results;
+        for (auto i : Range(curState.numEntries())) {
+            results.emplace_back(changes[i], make_pair(i, i));
+        }
+
+        const auto maxValue = max_element(results.begin(), results.end())->first;
+
+        return new ActivityAnimation(results, prevPositions.data(), curPositions.data(),
+                                     [=](float i) -> ActivityAnimation::Color {
+                                                 if (maxValue == 0) {
+                                                     return {1, 1, 1};
+                                                 } else {
+                                                     return {1 - i / maxValue, 1 - i / maxValue, 1};
+                                                 }
+                                             });
+    } else {
+        return new ImageInteractionAnimation(changes.data(), prevState.shape(), prevPositions, curPositions);
+    }
+}
+
+static Animation *getNormalizingAnimation(const fmri::LayerData &prevState, const LayerData &curState,
+                                          const vector<float> &prevPositions,
+                                          const vector<float> &curPositions) {
+    CHECK(prevState.shape() == curState.shape()) << "Shapes should be of equal size" << endl;
+    vector<DType> scaling(std::accumulate(prevState.shape().begin(), prevState.shape().end(), 1u, multiplies<void>()));
+    caffe::caffe_div(scaling.size(), prevState.data(), curState.data(), scaling.data());
+
+    // Fix divisions by zero. For those cases, pick 1 since it doesn't matter anyway.
+    normalize(scaling.begin(), scaling.end());
+
+    if (prevState.shape().size() == 2) {
+        EntryList entries;
+        entries.reserve(scaling.size());
+        for (auto i : Range(scaling.size())) {
+            entries.emplace_back(scaling[i], make_pair(i, i));
+        }
+
+        auto max_val = *max_element(scaling.begin(), scaling.end());
+
+        return new ActivityAnimation(entries, prevPositions.data(), curPositions.data(),
+                                     [=](float i) -> ActivityAnimation::Color {
+                                         auto intensity = clamp((i - 1) / (max_val - 1), 0.f, 1.f);
+                                         return {
+                                                 1 - intensity,
+                                                 1,
+                                                 1
+                                         };
+                                     });
+    } else {
+        transform(scaling.begin(), scaling.end(), scaling.begin(), [](float x) { return log(x); });
+        return new ImageInteractionAnimation(scaling.data(), prevState.shape(), prevPositions, curPositions);
+    }
+}
+
+Animation * fmri::getActivityAnimation(const fmri::LayerData &prevState, const fmri::LayerData &curState,
+                                       const fmri::LayerInfo &layer, const vector<float> &prevPositions,
+                                       const vector<float> &curPositions)
+{
+    if (prevPositions.empty() || curPositions.empty()) {
+        // Not all positions known, no visualisation possible.
+        return nullptr;
+    }
+
+
+    switch (layer.type()) {
+        case LayerInfo::Type::InnerProduct:
+            return getFullyConnectedAnimation(prevState, layer,
+                                              prevPositions, curPositions);
+
+        case LayerInfo::Type::DropOut:
+            return getDropOutAnimation(prevState, curState, prevPositions, curPositions);
+
+        case LayerInfo::Type::ReLU:
+            return getReLUAnimation(prevState, curState, prevPositions, curPositions);
+
+        case LayerInfo::Type::Pooling:
+            return new PoolingLayerAnimation(prevState, curState, prevPositions, curPositions);
+
+        case LayerInfo::Type::LRN:
+            return getNormalizingAnimation(prevState, curState, prevPositions, curPositions);
+
+        default:
+            return nullptr;
+    }
+}

src/fmri/visualisations.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include "LayerVisualisation.hpp"
+#include "LayerData.hpp"
+#include "Animation.hpp"
+#include "LayerInfo.hpp"
+
+namespace fmri {
+    /**
+     * Generate a static visualisation of a layer state.
+     *
+     * @param data
+     * @return A (possibly empty) visualisation. The caller is responsible for deallocating.
+     */
+    fmri::LayerVisualisation *getVisualisationForLayer(const fmri::LayerData &data, const fmri::LayerInfo &info);
+
+    Animation * getActivityAnimation(const fmri::LayerData &prevState, const fmri::LayerData &curState,
+                                     const fmri::LayerInfo &layer, const vector<float> &prevPositions,
+                                     const vector<float> &curPositions);
+}